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1. INTRODUCTION 


The character of natural buoyant convection in rigidly contained inhomogeneous fluids can be 
drastically altered by vibrating the container. Vibrations are expected to play a crucial influence 
on heat and mass transfer onboard the ISS. It is becoming evident (especially as demonstrated 
during the workshop held in ESTEC, September 98 [1]) that substantial vibrations will exist on 
the ISS in the wide frequency spectrum. 

In general, vibrational flows are very complex and governed by many parameters. It is 
almost impossible to correctly predict vibrational effects empirically; a sophisticated theoretical 
approach and numerical modeling are therefore essential. Available flight experiment data clearly 
show that, once initiated by “g-jitter”, convective flows can exist for a long period of time in a 
follow-up period of low gravity environment [ 2 - 1 ]. 

In many terrestrial crystal growth situations, convective transport of heat and constituent 
components is dominated by buoyancy driven convection arising from compositional and thermal 
gradients. Thus, it may be concluded that vibro-convective flow can potentially be used to 
influence and even control transport in some crystal growth situations. 

Control of convective transport continues to be an important aspect of crystal growth 
research. A control of convection through static and rotating magnetic fields has been used by 
several groups. In some cases, experimenters seeking to avoid buoyancy effects through use ot 
microgravity environments have expressed interest in the use of magnetic fields even under low 
gravity conditions. However, there are many instances, whether due to materials properties or 
other practical considerations, where the use of magnetic fields to induce stirring or suppress 
flow may not be an option. In such cases, vibrational control becomes an attractive alternative. 

Experimental results clearly show that vibrational convection can provide enhanced 
nutrient fluxes[8-14]. Numerical modeling work at the Institute de Mecaniques des Fluides de 
Marseille has confirmed that vibration can also be used to suppress buoyancy-driven flows. 
Furthermore, it was shown that such suppression would be very effective at reduced gravity 
levels of 10~ 4 g or less. This raises the exciting possibility that, under microgravity conditions, 
specific controlled vibration can be used to mask undesirable “g-jitter” induced convective 
effects. Such g-jitter convection can be caused by quasi-steady residual acceleration [15,16] (due 
to gravity gradient and atmospheric drag effects [7]), as well as transient and oscillatory 
acceleration disturbances [18,19]. It is understood that high-quality low gravity environments can 
only be provided for specific limited time intervals. For a space station, interruptions and 
disturbances are inevitable consequences of docking, pointing maneuvers, astronaut activity, etc., 
and will limit the maximum attainable duration of high quality microgravity periods. While 
active vibration isolation can be a partial solution, it will not solve the problems that might arise 
due to the quasi-steady and very low-frequency acceleration components related to the gravity 
gradient and other orbital factors. Alternatively, rather than using vibration to suppress low- 
gravity flows, one might envisage using vibration to provide flow regimes tailored to particular 
crystal growth experiments. These flows would not be accessible under terrestrial conditions due 
to strong natural convection effects. 

Bouyancy driven vibro-convective motion occurs when oscillatory displacement of a 
container wall induces the acceleration relative to the inner fluid. The vibration may be viewed as 
a time -dependent modulation of steady gravity in some cases. 
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In a closed container filled with a homogeneous fluid, the fluid moves as a rigid body 
with a container. If, however, the fluid density is nonuniform, fluid motion may ensue. Of course, 
the magnitude of this motion, depends on the orientation of the vibrational direction with respect 
to the local density gradients. Note that similar to Rayleigh-Bemard configurations there may be 
a “critical” threshold for the coupled vibrational frequency and amplitude, to cause convection. 
Interestingly, in the case of a constant density fluid, the fluid motion may also take place at 
angular vibration, axial rotational vibration or ACRT flows (caused by accelerated - decelerated 
rotation of a crucible) due to inertia forces. For density gradients caused by temperature, such 
motion is called thermo-vibrational. 

It has been recognized for some time that oscillatory or pulsatile flow can significantly 
alter the transfer of mass, heat and momentum in fluid systems [2,8-14,20-33], For example, 
analyses of heat transfer in laminar oscillating flows have shown that at high frequencies the 
effective diffusivity, k eff , behaves like k c({ ~ Ax 2 (tov) l/2 /L and, at low frequencies, k e{{ ~ 
Ajc 2 ( 0 3/2 v 1/2 /L, where v is the kinematic viscosity, Ax is the cross stream average of a fluid 
element over half the period of the oscillation and L is a characteristic geometric distance (e.g., 
between the plates). It was also shown that heat transfer was most enhanced when the 
characteristic heat transfer time was equal to half the oscillation period [24], 

A great deal of work related to the theory of thermo-vibrational convection has been 
carried out by Russian research groups. The main focus has been on thermo-vibrational 
convection starting with the work of Zenkovskaya and Simonenko [34], who first obtained the 
equations of thermo-vibrational convection in a high frequency limit. Since then, there have 
been many theoretical [2,35-42] and some experimental (e.g., [12,13,43]) studies of the stability 
of thermo-vibrational flows. One of the main conclusions that can be drawn from these works is 
that for vibrations with specially chosen axes, the natural buoyancy driven convection which 
would prevail in the absence of vibration can be suppressed at certain frequencies and amplitudes 
[39]. This has recently been analyzed in more detail [44] using the full equations of motion and 
Gershuni’s time-averaged equations. The possibility of using vibration as a means of controlling 
and suppressing convection was confirmed. A comprehensive introduction to vibrational 
convection and other time-dependent modulation can be found in reference [45]. 

The problem of vibrational convection arising due to other buoyancy sources, such as 
compositional density gradients, has also been approached. The onset of purely solutal and 
thermosolutal convection has been examined for horizontally stratified layers subject to vertical 
vibration [46,47], 

There are several examples of experimental work concerning the influence of vibration 
on crystal growth from melts and solutions [28-37]. These works involved a wide range of 
intensities and frequencies (including ultrasound). Experimental attempts to understand how 
low frequency vibrational stirring might be used to effect rapid mixing in melts and solutions 
have been made by Liu et al. [54], The influence of low frequency vibration on interface location 
and shape during Bridgman growth of cadmium telluride was examined by Lu et al. [55]. Other 
effects of low frequency vibrational convection on crystal growth include the increase in local 
perfection of binary compound semiconductors [57], changes in interface shape [30], and the 
facetting of germanium crystals [58], The elimination of striations in indium antimonide may 
also be due to the formation of a stationary melt flow due the torsional vibration [59]. 

Experimental results also clearly show that in certain cases vibrational convection can 
provide enhanced nutrient fluxes during the solution growth of Rochelle salt and potassium 
dihydrogen phosphate (KDP) [60,61], Zharikov [62, 14] identified a characteristic low 
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frequency (< 100 Hz) vibrational flow regime in the liquid near a growing crystal. The form, 
dimensionality and intensity of the flow were studied and the effects of vibration on heat and 
mass transfer were analyzed for the case of Czochralski and Bridgman growth of sodium nitrate 
(NaNC>3). He showed that the vibration could drastically alter the character of flow and 
concluded that vibration could exert a strong influence on transport and impurity incorporation 
and locally influence growth kinetics. 

Uspenskii and Favier [63] considered the interaction between high frequency and natural 
convection in Bridgman-type crystal growth. They used the average thermo-vibrational How 
equations to theoretically examined the problem of suppressing natural convection using high 
frequency (~ 10^ Hz) low amplitude vibration and compared the efficiency of vibrational 
damping to that of magnetic field damping. Using the physical properties representative of 
GaSb, GaAS, etc., they found that under terrestrial conditions, (for high electroconductivities) the 
magnetic field is more effective than vibration in damping flow in the horizontal Bridgman 
configuration. In contrast, for vertical Bridgman, lateral vibration was most effective. The 
horizontal velocity decreased by a factor of 10 and the vertical velocity by about 20. In 
comparison a 1 Tessla vertical magnetic field only resulted in a factor of 6 decrease in maximum 
velocity. They speculated that it might be possible to combine magnetic fields with vibration to 
reach optimal damping conditions. 

Vibrational flows are very complex and governed by many parameters. Here is an 
incomplete list: Pr, Ra, Rct v , Ra c , a r , co, UR (Prandtl, Rayleigh, vibrational Rayleigh, 

compositional Rayleigh numbers, a, = Ra v / Ra vr , and aspect ratio), q> , a (angles between 
gravity vector, vibrational direction and the axis (temperature gradient)), E, R 0 /L (Ekman 
number, relative distance to the center of rotation for angular, rotational vibrations, a type of 
vibrational action) and others. The range of parameter values is very wide: Pr is of the order of 
0.01 for semiconductors, and of the order of 10 for oxide melts; gravity variation result in Ra 
values of 10 6 to 10 10 on a ground to 1- 10 ? in space, and similar for Ra v ; the frequency range is 
from 0.0 1Hz to 100Hz and higher, etc. (see Table 1). 

Farooq and Homsy [41] examined in depth the 2D case of a rectangular cavity, but an 
understanding of 3D thermo-vibrational flows has not yet developed. Therefore, a sophisticated 
theoretical approach and numerical modeling together with experimental research are essential 
for the investigation of vibrational control of convective flows. 
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Table 1. Typical material properties for semiconductor and oxides melts, and 
nondimensional problem parameters for numerical experiments 

Dynamic parameters 


Frequency 
Angular amplitude 

Geometry 

Height 

Diameter 

Material properties 
Pr 

Viscosity 
Diffusivity 
Thermal diffusivity 
Thermal expansion coeff. 


0.001 ... 1 ... 100 Hz 
E = 0 ... 15 ...45° 


L = 2 to 10cm 

D = 1 to 2 cm => L/D= 10 to 5 


0.01 (Ge, Si) ... to 15 (oxides) 
v = 1.3e’ 2 cm^/s 
D = 1.3e" 4 cm^/s 
k = 1.3e 1 cm^/s 
P = 2.5e‘ 4 K-l 


Body forces and related parameters 


Gravity 

Temperature gradient 
Magnetic field 


g= 1 ...0.01 ... 10-5 ... 10-6 go 
G = 25 to 100 K/cm 

Ha = 0 ... 2,000 ... 20,000 


Nondimensional parameters 


v PATgL 3 P AT £2 0 L 4 ' v 

Pr=-,Ra = , Kan = , b = r, 

K KV KV QqL 

Ra v , L/D, h/R, Ha, Sc=v/D, Ra c *, ca= ^,/L 2 k 


Vibration types : 

translational, circular polarized, angular, g-jitter 


2. MODEL OVERVIEW 


2.1 Heat and Mass Transport in Bridgman melt configurations 

The research involved a numerical and experimental investigation of vibro-convective 
transport regimes applied to the control of convection and transport during crystal growth by the 
Bridgman technique. A typical Bridgman crystal growth configuration is shown in Fig. 1. This 
work was motivated by recent developments in the understanding of low frequency thermo- 
vibrational convection, results from recent Space Flight experiments on “g-jitter” induced flows, 
and by current and planned experimental work at Stanford University and the General Physics 
Institute, Russian Academy of Sciences in Moscow [1-5,40,64,65], Although in principle, the 
theoretical research described below can be carried out independently of a particular 
experimental program, close collaboration with experimental groups created a firm practical 
foundation. 

We investigated five types of vibro-convective flows caused by translational, angular, 
rotational vibrations (these results may be extended also for ACRT flows [92,94], due to 
accelerated-decelerated rotation of a crucible-tank), and flows due to “g-jitter”. The flow regimes 
were investigated by applying each type of vibrations to different materials to find regimes 
beneficial to control or suppress of convection. 

In section 2.2, the physical model and basic equations are described. The equations for 
translational, polarized vibrations and “g-jitter” are given in section 2.2.2. Section 2.2.3 describes 
the equations for angular, rotational and ACRT type vibrations. Solution methods are discussed 
in section 2.3. 

The basic problem analyzed was the suppression or control of buoyancy-driven 
convection in melts during plane-front directional solidification. The work involved the analysis 
of vibrational interaction with natural convection and its effects on the temporal evolution of 
melt temperature and composition during growth (section 3). The philosophy behind our 
approach is to use numerical modeling in two ways: synergistically, with experimental 
developments, and as a predictive tool. The synergism with experiment allows careful 
interpretation of both experimental and modeling results for what can be a highly non-linear 
physical situation. To explore regimes and system properties currently inaccessible to the 
experiment, the model was also used as a predictive tool. Initially guided by previous results 
[34-67,70,73-76], the work involved an extensive investigation of vibrational flow regimes with 
and without the presence of natural buoyancy-driven convection. As we progressed toward our 
goal of defining the applicability of vibrational control of convection, we compared the ability of 
vibration to suppress flow with that of a magnetic field. 
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2.2. Theory and numerical modeling 
2.2.1 The physical model 


The plane-front directional solidification of a binary solid from a melt of initial composition c„ 
in a cylindrical ampoule of radius R was considered. Fig. 1 . shows a physical model. 

Calculations were made for two basic physical set-ups. 

The first involves purely 3D thermo-vibrational convection in a differentially cylindrical 
cavity with no consideration of solidification. These calculations were performed in conjunction 
with the flow visualization experiments, ongoing experiments conducted by R.S. Feigelson at 
Stanford University. 

The second setup is closer to a practical solidification situation and involves explicit 
consideration of plane-front solidification in a long cylindrical ampoule. Solidification takes 
place as the ampoule is translated along a temperature gradient. For this model system, 
translation of the ampoule is simulated by supplying a doped melt of bulk composition at a 
constant velocity V m at the top ot the computational space (inlet), and withdrawing a solid of 
composition c s = c s (x,t) from the bottom. The crystal-melt interface is located at a distance L 
from the inlet; the radius of the ampoule is R. The temperature at the interface is taken to be T m , 
the melting temperature of the crystal, while the upper boundary is held at a higher temperature 
7), (Fig. 2). On the ampoule walls the heat transfer between the furnace and the ampoule will be 


HEATER 


COOLER 


TRANSLATION 

DIRECTION 



MELT 


THERMAL 

BARRIERS 


CRYSTAL 

AMPOULE 


Fig. 1 The Bridgman-Stockbarger method in which 
the crystal is directionally solidified by 
translating an ampoule between hot and cold 
zones of a furnace. 



Crystal 

2R 

Fig. 2. The model Bridgman system. V m is the 
ampoule translation rate, T a (z) is the heater 
profile and T m is the melting temperature 
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modeled using a suitable heat transfer condition. The heating profile includes a long isothermal 
hot zone. In this case, thermo-vibrational convection in the upper portion of the melt will be 
weak or nonexistent. Thus, the effect of our implicit neglect of the no-slip condition at the top of 
a real ampoule will be negligible. 

Both thermal, solutal and thermo-solutal convection were considered. For thermo- 
vibrational convection we considered two types of vibration, translational and angular. The frame 
of reference for the calculations is an experimental frame attached to the walls of the ampoule. 
Thus, depending on the type of displacement induced by the vibrating walls, the equations 
governing momentum mass and heat transport will take on a different form. In section 2.2.2 we 
consider translational and polarized vibrations, followed by angular and rotational vibration in 
section 2.2. 3. 

2.2.2. Translational, polarized vibration and “g-jitter” 

Translational vibration corresponds to a linear displacement such as, for example, u = 
djeoscot , where d \ is a real vector corresponding to the displacement amplitude. Thus, a point is 
displaced back and forth upon the same line. Polarized vibrations are characterized by a 
displacement u — Re(de~^ w ^ } where d = dj - / d2 (see Fig. 3). Here the instantaneous vibration 
direction rotates in the polarization plane defined by the real vectors dj and r/2- A case of 
orthogonal equal length vectors dj and rf? , the circular polarized vibration, is called a CVS 
method [32] 



Figure 3: Translational vibration (a), d j or d2 = 0, and polarized vibration (b), d},d2 ^ 0; <p is the angle 
between gravity vector and axis of ampoule, OC As the angle between direction of vibrations and axis of ampoule 

For a reference frame fixed to the vibrating ampoule both these types of vibrations result 
in equations of the following form [17,46,47]. The momentum and continuity equations are 

A + (y.V)V + Vp-PrV 2 V = RaPr(Q + ac)k + Ra v Pr(Q + a r c)f(\,0),t) (1) 

a 

VV = 0 (2) 

while the continuity and heat - mass transport equations are: 
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( 3 ) 


— + V • V0 = V 2 0 
dt 

■j 

— + V • Vc = Pr Sf-'v 2 c (4) 

dt 

where length, time and velocity have been scaled by Rq, R(^/k and k/Rq, where Rq is the radius 
of the ampoule and K is the thermal diffusivity, k is a unit vector in the direction of gravity. 


The Prandtl, Schmidt, Rayleigh and solutal -thermal buoyancy ratio are given by 



Sc= Ra = 


P* TgR 3 

VK 


_p c c„ 

^ “ PA T ’ 


where P and p s are the thermal and solutal expansion coefficients and AT, g, v and D are the 
characteristic longitudinal temperature difference in the melt, gravitational acceleration, 
kinematic viscosity and solute diffusivity, respectively. The dimensionless number Ra v and Ra vc 
are vibrational Rayleigh numbers defined as 


_ d w 2 PA tr> 


^ 2 o A ~r > 3 


For the solidification calculations, the following boundary conditions apply at the crystal-melt 
interface. 


0=QV-N=P<?,VxN=O 


(5) 


— = P eSc Pr~'(l - k)c 
dn 

wherei’e = V m R/k is the Peclet number, k is the distribution coefficient and N is the unit vector 
normal to the planar crystal melt interface. At the top of the ampoule, 


^=l,V-N=Pe,VxN=0 


( 6 ) 


— = PeSc Pr"'(c - 1) 

dn 

where, due to the fact that we consider the top of the computational space as an inlet, the 
last condition represents conservation of mass. Finally, on the ampoule walls we have 

V-N = Pe g Sc Pr~‘, V-e w = 0, We w = h(0-^(z)), VC-e w = 0, (7) 

where © a (z) is the heating profile, h is a dimensionless heat transfer coefficient and e w is 
normal to the walls. 
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The body force f(x,r) = nsin(£20for translational vibration, where n is direction of 
vibration (Fig. 3b); f(x,t) = [cos( Q. t), sin( Q t), 0] T for polarized vibrations and f(x,t) = J g (t), the 
3D vector of “g-jitter” microaccelerations (e.g. measured SAMS microacceleration data or the 
numerical “g-jitter” model). 

Vibrating immersed body. The application of vibration through placement of vibrating 
bodies at specific locations in the melt can also be modeled as a variant of the translational and 
polarized vibration problems. 

Average equations for small amplitude vibrational displacements 

The dimensionless groups Ra v and Ra vc in equation (1) should not be confused with the 
vibrational Rayleigh number Rcig of Gershuni [34-39,45-76], which arises for monochromatic 
vibration as follows. If the amplitude |d| of a sinusoidal vibration is sufficiently small, i.e., |d| « 
min( R/PAT, R/ Pcoo), the velocity, temperature and concentration fields may be represented as 
the superposition of mean (averaged) fields and small oscillating components [34,74], In this 
case, application of the averaging method [45] yields 

^ + (rV) v = - Vp + Pr V 2 r + PrRa( 0+ c ^ C) k 4 -PrRa g (»wV)[(0+ C) n - w], (8) 

V • v = 0 , (9) 

where Rag = (Pen d RAT)2 /kv is Gershuni’s thermo-vibrational Rayleigh number and = 
PcCoo/PAT is the ratio of the solutal and thermal buoyancy. The vector w is the slowly time- 
varying part (may be stationary) of the oscillating part of the velocity field and satisfies 

V • w = 0 , Va w = V[(0 + a c C)] An , (10 ) 

and is the solenoidal part of the vector field (0+ a c C)n = w + V<p . The average temperature 
and solute fields satisfy 

+ (r-V ) 0 = V 2 © , ( 11 ) 

and ' _ 

+ ( f V ) c = |t-r ! c 02) 

The boundary conditions for the averaged equations are of the same form as (4)- (8) 
and the component of w normal to the boundaries must be equal to 0 on the ampoule walls and 
equal to Pe at the top of the ampoule and at the melt-crystal interface. For thermally buoyant 
natural and vibrational convection we simply set etc = 0 and for no solidification, Pe = 0. 

For the problems under consideration, the averaged equations result in the steady-state 
equations. Our recently developed novel methods of numerical bifurcation theory [69] are 
suitable for the analysis of stability and onset of vibrational flow regimes. 
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Figure 4: Diagram of the system for angular vibration to be modeled (a), x* and x are the two coordinate systems 
related by (10). The container is rotated at an angle T(t) = f sin ( Q 0 t) about a center of rotation at x* = 0. 
Axial vibrations and ACRT (b), are rotation about the axis z, the angle T(t) = f sin ( £l 0 t) for vibrations or 
T=T(t) (selected ACRT law). 



2.2.3. Angular, rotational vibration and ACRT flows 

The equations of motion for angular and rotational vibrations take on a more complicated form. 

Angular vibration. Let x* = (x [*, X2*, X3*) be a fixed Cartesian coordinate system (see Fig. 4a). A 
container of length L is subjected to an angular vibration q(t) in the xj*-x 3* plane. The position 
vector to the mass center of the cylinder is parallel to the side of the cylinder and is given by 

q* = Rq[- sin 0q + cos Brj], (13) 

where Rq is the distance from the origin 0 * to the mass center of the cylinder and 

0(t) = esin Q 0 t , (14) 


where t is time. 

Rotational vibration and ACRT (Fig. 4b). The container is rotated about the axis z, the angle T(t) = e 
sin (Q. 0 t) for vibrations or T=T(t) (selected ACRT law, see [77]); <p is the angle between the gravity 
vector and axis of ampoule. The ACRT type of vibration was not in the original research plan. We 
just note here, that the software developed in the frame of current research is capable for the 
investigation of ACRT flows without any modifications. 

In the reference frame moving (rotating) with the ampoule the momentum equation has a form 


— + (V • V)V + p~ x Vp - vW 2 \ = yTg + (CO x q)x a> + 2V x co + q x rff (15) 
dt 

where co= [0, e d 0 cos(Q. 0 t), 0] T for angular vibration and (O = [0, 0 , e Q 0 coj(£2 0 f)] T for rotational 
vibration and co= [0, 0, dT(t)/dt ] T for ACRT. The problem has the following nondimensional 
determining parameters: Q = QoL 2 /k, the dimensionless frequency and 0 = Ro/L; the Prandtl, Pr, 
Rayleigh, Ra, vibrational Rayleigh Ra q, and Ekman, E , numbers are 


Pr = Ra 


P AT gL 3 

KV 


Ra 


a 


P AT Q0L 4 

KV 


, E = 


v 

£2 0 L 2 ’ 


where P, v, g and K are the coefficient of thermal expansion, kinematic viscosity, gravitational 
acceleration and thermal diffusivity, respectively. Note that, the system of equations (15) differ from 
the usual equations in the absence of rotation in that two additional terms are present in each 
equation; the Coriolis term which is proportional to zPr/E , and the centrifugal term which is 
proportional to E^ORa^Pr at x = 0 and varies with position in the ampoule. The importance of the 
latter term depends on the dimensions of the amplitude of the angular vibration, E , and the ratio 0. 
Note that the centrifugal terms give rise to periodic forcing that fluctuates about a mean value at twice 
the period of the angular vibration. 

Since the above system of equations has not been well studied, we used a conservative 
approach for the study of angular vibrations and confineed our investigation to a parametric study of 
flow regimes and transitions for thermo-vibrational situations only, without solidification. 
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2.3. Solution method and software 

The equations are solved in primitive variable form (velocity-pressure, temperature, etc.) using a 
Finite Element Method package developed by the Principle Investigator [65-67,71-73]. The package, 
FEMINA/3D , is designed to solve 2D and 3D problems and includes four basic groups of modules 
each distinguished by its function: automatic 2D/3D mesh generation, optimal renumbering or 
matrix-bandwidth/front optimization, a computational part including efficient algebraic solvers and 
visualization/graphics software. The novel methods developed and implemented allow a use of 
efficient, same order finite element approximation, for velocity and pressure, and robust iterative 

solution techniques of a high accuracy (up to 10 ’). 

The FEMINA/3D code has overcome most of the disadvantages associated with the memory 
requirements of FEM methods by using an efficient combination of solution methods. These include 
the iterative solution of large non-symmetric algebraic systems using iterative IDR-CGS type 
techniques [72,82-84,91] with preconditioning by incomplete decomposition of high order in a 
compact matrix storage scheme. That reduced the computation time by one to two orders of 
magnitude and memory requirements by an order for 3D flows over currently available commercial 
codes (e.g. CFD2000). Typical solution time for a transient 3D problem with a 200,000 unknowns is 
a couple of hours on a SGI 02 workstation (much faster on CRAY-C90); solution time for a 3D 
steady convective flow problem is a couple of minutes. 

FEMINA/3D simultaneously considers the continuity equation directly with the momentum 
equations at each solution step. This eliminates many problems related to boundary conditions and 
places only slight limitations on the time step size for transient problems (these limitations are related 
to the physical nature of the problem rather than the numerical method itself). 

The package has been carefully verified by comparison of numerical and analytical solutions 
as well as comparison with results obtained by other methods and experimental data for 2D and 3D 
viscous flows and thermal convection problems [65,67,73,86,91]. 

A wide range of problems has been investigated using this package, including time-dependent 
fluid flow (64, see also [72] and references therein) and mass and heat transfer in complex geometries 
[64,78,79], Other practical applications include thermal convection and transport during crystal 
growth processes [33,40,79,80,81] and other engineering problems. 
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3. INVESTIGATION OF VIBRATIONAL CONTROL OF CONVECTIVE FLOWS IN 

BRIDGMAN MELT CONFIGURATIONS 


Results of the investigation of vibrational control of the Bridgman crystal growth technique have been 
published in journal papers, conference proceedings and presented at conferences, workshops and 
colloquia, listed in Section 6. 

A summary of results from our publications [33,40,64-70,86-92,99] is presented below. 

(i) An efficient numerical code, a main research tool, for the modeling of three-dimensional 
thermal-vibrational convection for several types of vibrations has been developed by the PI based on 
the finite element code FEM1NA/3D created by the PI [71,72], The code was carefully tested on 
experimental benchmarks, and published numerical data for a variety of 2D and 3D viscous flows and 
thermal convection problems [65,67,73]. The novel methods developed and implemented allow a use 
of efficient, same order finite element approximation, for velocity and pressure, and robust iterative 

solution techniques of a high accuracy (up to 10~ 9 ). This reduced the computation time by one to two 
order of magnitude and memory (by 8 times less) lor 3D flows over currently available commercial 
codes (e.g. CFD2000). Typical solution time for a transient 3D problem with 200,000 unknowns is a 
couple of hours on a SGI 02 workstation; solution time for a 3D steady convective flow problem is a 
couple of minutes [86,90,91]. A new efficient numerical continuation approach, based on the 
multiquadric method, for the stability and bifurcation analysis was developed [69, 89]. 

(ii) A parametric study of the general 3D buoyant-vibrational flow in a Bridgman growth 
system has been performed. This includes the characterization of flows induced by translational, 
rotational and circular polarized (also known as CVS, after R. Feigelson [32]) vibrations under 
selected microgravity and Earth conditions for typical semiconductor and oxide melts [33,40], A 
typical flow pattern for translational vibration is presented in Fig. 5. Even in Og we have found flows, 
i.e. characterized by a fluid flowing up one side and down the other (relative to the Figure). The 
temperature distribution remain almost unperturbed for semiconductor melts (due to the low Pr 
number), therefore a separate plot for AT is provided. 



Vx Vy Va P T AT l|V|| 

Figure 5. Typical instantaneous 3D melt flow patterns for translational vibration at Og, Ra=0, Ra r =7 . 25 • 10 4 , (0 = 100Hz, 
lateral vibration: velocity components V„V y ,V z , pressure P, temperature T, temperature disturbance AT, and velocity 
module |V|. Black color designate the minimal value plotted, white one - the maximal value 
The angular orientation between the direction of vibration and ampoule axis (with imposed T- 
gradient) has been studied for translational vibrations. When that angle is zero, there is no influence 
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of the vibration on a flow even if vibrational Raleigh number is very high. The maximum observed 
effect corresponds to an angle of 90 degrees. 

The flow patterns for rotational vibrations flow regimes are presented in Fig. 6. Maximal 
values of flow velocity are observed at the end of the ampoule which is more far from the rotation 
origin. 



t A T uni 


Figure 6. Typical instantaneous 3D flow patterns for angular vibration at Og, Ra=0, Ra v = 4.6 • 10 , (0 - 100Hz, lateral 
vibration: velocity components V x ,V y ,V 2 , temperature T, temperature disturbance AT, and velocity module |V|. 

We found that both translational and angular vibration can cause average melt flow for a range of 
parameters typical of practical semiconductor and oxide growth. For a given vibration amplitude and 
frequency, angular and circular polarized vibrations result in more intensive melt flow than 
translational ones. 

(iii) The influence of vibrations on the heat - mass transfer becomes much more significant for 
oxide melts due to their low thermodiffusivity (Prandtl number ~ 10, that is four orders higher than 
for semiconductors). These flow patterns are shown in Fig. 7 for the case of circular polarized 
vibration. 



Vx Vz T C |V| A T 


Figure 7. Instantaneous 3D flow patterns for circular polarized vibration, Ra = 7.25 • 10\ Ra „ = 7.4 ■ Iff , 0) - 10Hz : 

velocity components V x , W z , temperature T, concentration C, velocity module |V|, and temperature disturbance AT. 
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Initially (at the time = 0 ) the species concentration was c-l at the lower quarter of the cylinder and c 
- 0 elsewhere. The evolution of the species concentration (process of mixing) is shown in Fig. 8 
together with minimum and maximum values of velocity components. Complete mixing 



Figure 8. Temporal evolution of concentration C (left, min. and max. of Cj and velocity extremums (V x , center, and V z , 
right) for circular polarized vibration, Ra=7.25 - 10\ Ra^-IA- 10 6 , (0 = 10Hz. 


occurs in about ten seconds. The heat transfer (local Nusselt number at the top and the bottom) is also 
enhanced by about an order of magnitude. If the frequency of vibration is higher, of the order of 
100Hz (for fixed Ra v ), then the change of heat and mass transfer due to vibrations becomes less 
significant. That corresponds to the experimental observations by R.Feigelson [32,74]. The 
visualization of results of numerical modeling, the 3D movies of vibrational flows, heat and mass 
transfer are available on the Web, http://uahtitan.uah.edu/alex/cvs_numeric.html for some cases. 

(iv) The oscillation of flow values, in addition to a drive frequency, shows co/2 subharmonic 
in the pressure field. Similar experimental spectra was observed by M.Schatz [75]. With the 
increasing of vibration intensity, additional oscillating flow frequencies (superharmonics) appear. 

(v) The results on the influence of forced vibration on g-jitter induced flows are available 
[64], G-jitter was implemented using SAMS microacceleration data from the USML-2 (Fig. 9). 
Initially the translational vibration were applied in the direction parallel to the ampoule axis (T- 
gradient), trying to damp flow variation with time caused by g-jitter (predicted by the averaged 
equation theory [45]). 

USML-2 DAY 1 HEAD A SAMS DATA 



Figure 8: G-jitter microacceleration SAMS data from USML-2 mission, used in numerical simulations. 
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While the flow variation with time becomes more regular, a suppressing temporal variation of the 
melt velocity due the g-jitter was not succeeded. Alternatively, the use of the same amplitude 
vibration in the direction orthogonal to the ampoule axis induces intensive thermal vibration flows 
and flow disturbances due to g-jitter become practically not visible. Nevertheless recent results show 
that axial translational vibrations have a tendency to prevent the development of a strong stationary 
vortex cell, that otherwise may occur. 

(v) Convective flows in a semiconductor melt with strong static axial magnetic field applied 
were investigated, and we compared three different numerical methods for the solution of thermal 
convection [87]. Although the generated flows are extremely low, the computational task is very 
complicated because of the thin boundary layer at high Hartmann numbers, Ha » 1. We considered 
melt region geometry with different aspect ratios, and gravity direction aligned and misaligned with 
the magnetic field vector. The comparison shows that the finite element approach with regularization 
can obtain stable and reliable solutions in a wide range of Ha number, up to 10 . These results 
compare favorably with asymptotic solutions [99,100]. Numerical solution of these problems by 
available commercial CFD codes may be not efficient or not possible at all. 



(a) 



THERMAL CONVECTION - VERTICAL VELOCITY 



(C) 



VERTICAL VELOCITY IN THE BOUNDARY LAYS Ft 



Figure 9:Stream function for aligned (a) and (b) misaligned by 0.5 degree gravity and magnetic field directions, Ha=2170 
(B = 5T1) (c) Vertical velocity profile Vy(x) at y=0.25H, x in [-1,1], (note boundary layer), Ha=2.2 10 4 , B = 50 T1 (d) 
Detail of the velocity profile (c) in a boundary layer, x from -1 to -0.99. (e) Summary of the magnetic field 
suppression of melt flow: maximum value of horizontal velocity magnitude versus magnetic field B for aligned and 
misaligned by 0.5 degreemagnetic field and gravity vectors. Predicted theoretical asymptotic dependence for the 
velocity Vm^-Ha 2 is observed from about B=0.5 Tl. Results are obtained with our FEMINA/3D code for Ha number 
up to Ha=2 10 4 ; other tested in methods did not provide acceptable results or failed for Ha > 100 
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The efficiency of the convective flow suppression by axial magnetic field at Ig is very sensitive 
to misalignment of the gravity vector with ampoule axis direction. For example, in case of 
misalignment by 0.5 deg., the axial temperature gradient (which was stabilizing) becomes a main 
driving force for the thermal convection, resulting in the increase of convection velocity 
magnitude by a factor of 3 to 10. 


4. CONCLUSIONS 

• We verified and confirmed the validity of the Boussinesq model for semiconductor and 
oxide melts under microgravity conditions for all the flow regimes under consideration. 

• Powerful tool, the 3D finite element code FEDIMA/3D has been developed. Efficient 
numerical methods for 3D thermal vibrational convection are implemented; carefully tested; 
good agreement with experimental data is obtained. This tool is capable tor high accuracy 
modeling of 3D thermo-vibrational convection, and flows under strong magnetic field in a wide 
range of problem parameters. 

• A parametric study of the general buoyant-vibrational flow in a Bridgman growth system 
was performed for both semiconductor and oxide melts. 

• The influence of angle between a direction of vibration and ampoule axis (temperature 
gradient) has been studied for translational vibrations. Zero angle - no influence of the vibration 
on a flow. The maximum effect corresponds to an angle of 90 degrees. Here transport is 
significantly enhanced. 

• Influence of g-jitter and forced vibrations has been analyzed. 

• All kinds of vibration can cause average melt flow for a parameter range for praclical 
semiconductor/oxide growth. Rotational (angular) and polarized vibrations result in more 
intensive melt flows than translational ones. Typical flow patterns for different flow regimes 
have been identified. 

• The vibrational influence is stronger for oxides than for semiconductors. 

• The frequency range was identified where vibrational influence efficient. 

• Suppression of convective flows in a semiconductor melt with strong static axial 
magnetic field applied were investigated High sensitivity of magnetic field direction and gravity 
vector misalignment has been found and analyzed. 
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Abstract 

It is .-enenllv recognized that oscillators', or pulsatile, flow significantly alters the transfer of mass, heat and 
momentum n fluid s/stems. A numerical investigation of thermovibrational buoyancy-driven flow in different a > 
— jmitc’l containers is presented as part of a study of thermovibrattonal transport regimes m Bndgman-type 
systems The formulation of a physical and mathematical model for this problem is outlined and its application to the 
study of investigation of thermal vibrational flows is discussed. Three types of v, brat, on are 

circularly polarized and rotational. It is demonstrated that forced vibration can .significantly ««*« ^ V SOM EhS 
induced by u -jitter and selected results for the cases of longitudinal and lateral vibrations are presented. ^ - 

Science B.Y. All rights reserved. 
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1. Introduction 

The character of natural buoyant convection in 
ricidly contained inhomogeneous fluids can be 
drastically altered by vibration of the container. 
For certain experiments and operating conditions, 
vibrations are expected to have a significant influ- 
ence on heat and mass transfer onboard the Inter- 
national Space Station (ISS). Furthermore, it 
appears that (see for example, the recent ESTEC 
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Workshop proceedings [1]) (/-jitter vibrations will 
exist on ISS over a wide range of frequencies. 

In general, vibrational flows are very complex 
and are governed by many parameters. This com- 
plexity makes it almost impossible to correctly pre- 
dict vibrational effects empirically. Thus, a careful 
theoretical approach combined with numerical 
modeling is essential. Available flight experiment 
data clearly show that, once initiated by “(/-jitter”, 
the effects of convective flows can persist for long 
times even when the (/-jitter disturbance (and con- 
sequent flow) were short-lived [2-7], 

In many terrestrial crystal growth situations 
convective transport of heat and constituent com- 
ponents is dominated by buoyancy-driven convec- 
tion. Control of convective transport continues tc 
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be an important aspect of crystal growth research. 
Several groups are actively pursuing control of 
convection using static and rotating magnetic 
fields. In some cases, experimenters seeking to 
avoid buoyancy effects by conducting their experi- 
ments in microgravity environments have ex- 
pressed interest in the use of magnetic fields even 
under low gravity conditions. However, there are 
many instances, whether due to materials proper- 
ties or other practical considerations, use of mag- 
netic fields to induce stirring or to suppress 
unwanted flows may not be an option. Magnetic 
fields cannot be used for flow control in melts and 
solutions that are poor conductors. Flow sup- 
pression through vibration or vibro-convective 
mixins mav offer an attractive alternative in such 
cases. 

It is understood that high-quality low-gravity 
environments can only be provided for specific lim- 
ited time intervals. For space stations, interruptions 
and disturbances are inevitable consequences of 
docking, pointing maneuvers, astronaut activity, 
etc., and will limit the maximum attainable dura- 
tion of high-quality microgravity periods. While 
active vibration isolation can be a partial solution, 
it will not solve the problems that might arise due 
to the quasi-steady and very low-frequency acceler- 
ation components related to the gravity gradient 
and other orbital factors. Alternatively, rather than 
using vibration to suppress low-gravity flows, one 
micht envisage using vibration to provide flow re- 
gimes tailored to particular crystal growth expeii- 
ments. 

Recent work has shown that the character of 
natural buoyant convection in nonuniformly 
heated, rigidly contained inhomogeneous fluids can 
be drastically altered by vibration of the container. 
A review and relevant theoretical and experimental 
research can be found in publications [1-13]. Thus, 
vibrational, induced flow can potentially be used to 
influence and even control transport in some crys- 
tal growth situations. A practical quantitative un- 
derstanding of vibrational convection as a control 
parameter in crystal growth situations is currently 
not available. The objective of the work is to assess 
the feasibility of the use of vibration to suppress, or 
control, convection in order to achieve transport 
control during crystal growth. 


2. Theory and numerical modeling 

Buoyancy driven vibro-convective motion oc- 
curs when oscillatory displacement of a containei 
wall induces the acceleration of a container wall 
relative to the inner fluid. The vibration may be 
viewed as a time-dependent modulation of steady 
gravity. In a closed container the fluid will move as 
a. rigid body with a container. If, however, the fluid 
density is nonuniform, fluid motion may ensue. The 
magnitude of this motion, of course, depends on the 
orientation of the vibrational direction with respect 
to the local density gradients. Note that, similar to 
Rayleigh-Benard configurations, there may be 
a “critical” threshold for the coupled vibrational 
frequency and amplitude, to cause convection. 
Interestingly, it should be noted that in case of a 
constant density fluid subjected to spatially non- 
uniform vibration, fluid motion can also occui (foi 
example, angular vibration [11]). 

To properly investigate influence of transla- 
tional. circularly polarized and rotational (angular) 
vibration necessitates the use of the full 3D equa- 
tions governing the transport of heat, mass and 
momentum. Selected examples of our ongoing 
work on this topic are outlined below. 

2.1. The physical model 

We consider a purely thermo-vibrational con- 
vection in a differentially heated cylindrical cavity 
with no consideration of solidification. The fluid is 
taken to be Newtonian, with a constant viscosity 
and the Boussinesq approximation is assumed to 
hold. The validity of this approximation is dis- 
cussed in Section 3. The calculations were per- 
formed for identification and characterization o: 
thermovibrational flow and are part of an ongoing 
project involving flow visualization model experi 
ments being conducted by Feigelson [10]. 

2.2. Mathematical model for translational vibratwi. 

Translational vibration corresponds to a linea 
displacement such as. for example, u = d cos co, 
where d is a real vector giving the displacemer 
magnitude and to is the frequency. In this case th 
ampoule is displaced back and forth upon the sam 
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Fig, 1. Translational vibration (a). d x 
ampoule axis, u is the angle between 


or d, = 0. and polarized vibration (b). d,.d 2 is the angle between gravity vector and the 

the vibration direction and the ampoule axis. 


line. Polarized vibrations are characterized by 
a displacement it = Re{</e ia>r } where d — d { — i d 2 
(see Fie. 1). Here the instantaneous vibration direc- 
tion rotates in the polarization plane defined b% the 
real vectors d 1 and d 2 . A sketch showing both 
translational and circular polarized vibrations is 
presented in Fie. 1. In a refeience fiame fix^d to 
a vibrating ampoule, these types of vibrations result 
in the following form of the momentum equation: 

— + (r • P)r = -Vp + Pr V 2 v + Ra T Pr(0 + jClft 
5 t 


4- RafPr(0 4- y.C)f{Q.i). (11 

where length, time and velocity are scaled by 
R 0 , Rz/k and k/R 0 . Here R 0 is the ampoule radius 
and k is the thermal diffusivity. The nondimen- 
sional concentration and temperature, are given by 
0. and C. respectively. The function f{Q. r ) is the 
acceleration of the vibrating ampoule and 
Q = coRI/K is a dimensionless frequency. The con- 
tinuity and heat-mass transfer equations complete 
the problem formulation. The Prandtl. Schmidt, 
thermal and solutal Rayleigh and vibrational 


Rayleigh numbers and the buoyancy ratio are. re- 
spectively, given by 


Pr = — . Sc = — . Ra T 
k L> 


VK 


/I C t TO 

= PAT' 


Ra* 


dor/jAPR 3 

VK 


(21 


Here /I and j? c are the thermal and solutal expan- 
sion coefficients and AT. c x , g, d. a>, k, v, D are the 
characteristic longitudinal temperature diffeience. 
reference concentration in the melt, gravitational 
acceleration, vibrational displacement amplitude 
and frequency, direction of gravity, kinematic vis- 
cosity and solute diffusivity, respectively. The di- 
mensionless number Ra^ is the vibiationa 
Rayleigh number and Raf = aRa*- Eq. (21 is solvec 
together with the equations governing heat anc 
species transfer and the condition that the velocity 
is divergence free. 


2.3. Rotational vibration 

The equations of motion for angular vibration 
take on a more complicated form (see Fig. 2 
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Fie. 2. Diagram of the rotational (angular) vibration. Tr.e con- 
tainer is rotated at an angle <)(r) about a center of rotation at 
v * _ o jh e vector q* connects the center of rotation to tr.e mass 
center of the container. 


A container of length E is subjected to an angulai 
displacement 0(f) in the plane. Heie the cooi- 

dinates at* are referred to a fixed laboiatoiy trame 
of reference. The position vector to the mass ^entei 
of the cylinder is parallel to the side of the cylinder 
and is given by q* = RoL ~~ s * n ^*i + cos 
where R 0 is the distance from the origin 0* to the 
mass center of the cylinder and 0(r) = e sin Q 0 t. In 
a frame of reference moving with the container, the 
equations of motion have the form 



div v — pA'(sin Qi\ + cos 0/ 3 ) 


+ p\_2i2v - Q'x + Q.x 


-t- £~£2qRo COS 


(3) 


where D/Dr denotes the material derivative, r is the 
velocity of the fluid relative to the moving reference 
frame, p is the density of the fluid. £2 is the .ate of 
rotation tensor for the moving frame with respect 
to the fixed frame of reference. £2 is its time deriva- 
tive and T is the Newtonian stress tensor for the 
fluid. As in the previous example, the fluid is taken 
to have a constant viscosity and the Boussinesq 
approximation is assumed to hold. The dimension- 
less equations governing the transport o! mo- 
mentum. mass and heat in the cylinder are obtained 
after using L. L~/k. k/L, and AT = Th Tq to 
scale, respectively, length, time, velocity and tem- 
perature. The governing dimensionless parameteis 
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are the dimensionless frequency Q = QqL-’k. the 
dimensionless container radius, 5=K 0 /T. the 
Prandtl, Pr. and the Rayleigh. Ra. vibrational 
Rayleigh Ra a , and Ekman, E. numbers. The latter 
are given by 



„ 0A TgL> 

Ra = 

KV 



pATQlL* p 

KV ’ QoL 2 


(4) 


where p, v, g and k are the coefficient of thermal 
expansion, kinematic viscosity, gravitational accel- 
eration and thermal diffusivity, respectively. This 
system of equations differ from the usual equations 
in the absence of rotation in which additional teims 
are present; the Coriolis term which is proportional 
to ePr/E, and the centrifugal term which is propor- 
tional to e 2 5Ra fl Pr and varies linearly with position 
in the ampoule. The importance of the lattei term 
depends on the dimensions of the amplitude of the 
angular vibration, e, and the ratio 9. The rocking 
motion of the angular vibration under considera- 
tion means that centrifugal terms give rise to a peri- 
odic forcing that fluctuates about the mean value at 
twice the period of the angular vibration. 

Since the above system of equations has not been 
well studied, a conservative approach was adopted 
for the study of angular vibrations and we confine 
our investigation to a parametric study of flow 
reaimes and transitions for thermo-vibrational situ- 
ations in the absence of solidification. 


2.4. Solution method 

The equations are solved in primitive variable 
form (velocity-pressure, temperature, concentra- 
tion, etc.) using a finite element method package 
FEMINA/3D [14]. The continuity equation and 
momentum equations are consideied simulta- 
neously at each time step. This eliminates many 
problems related to boundary conditions and pla- 
ces only slight limitations on the time step size-foi 
transient problems (due to the physical natuie ol 
the problem). The regularization for the incompres- 
sibility condition makes the solution proceduit 
more efficient, and allows the same order finit< 
element approximation for both the velocity anc 
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pressure [15]. This approach makes it possible to 
solve large time-dependent problems (up to 300 000 
unknowns) on a SGI 02 workstation with reason- 
able computation times. 

We implemented the above 3D models of 
convective buoyancy-driven melt flow in differen- 
tially heated cylindrical containers using the 
FEMINA/3D code. This code was carefully tested 
on benchmarked experimental and numerical data 
for a variety of 2D/3D viscous and thermo-con\ec- 
tive flow problems [15.16]. 

For rotational vibrations the Ekman number 
can be of the order 10“ 4 -10" s for frequencies on 
the order of 1 Hz. This results in large coefficients, 
p r ■/£, for the Coriolis terms in the governing 
equations and causes difficulties in the numerical 
solution. To resolve this we implemented a high- 
accuracy solution method using preconditioning 
by high-order incomplete decomposition. This al- 
lowed us to obtain high-precision solutions with 
accuracy up to 10" 9 . Preconditioning also reduced 
the computation times by one to two orders of 
magnitude and the memory size by a factor of S for 
3D flows compared to currently available commer- 
cial codes (e.2. CFD2000). A typical solution time 


for a transient problem is about 2 h on a SGI 02 
workstation. 

3. Results and discussion 

We verified the validity of the Boussinesq model 
for semiconductor and oxide melts under micro- 
gravity conditions. This topic has been discussed 
recently by Perera and Sekerka [17], Pukhnachev 
[18] and Gershuni and Lyubimov [11]. If the non- 
dimensional criteria, proposed by Pukhnachev. 
Pu = gL 3 v _1 K~ 1 is less than 1, then the Boussinesq 
model for thermal convection may not be valid. 
Our estimates show that the Boussinesq model is 
quite adequate for the differentially heated closed 
ampoule and the range of parameters and material 
properties under investigation. The values of Pu are 
of the order 10 4 -10 5 for semiconductor and oxide 
melts for g/g 0 = 10 _5 -l(T 4 , clearly well above 1. 

A parametric study of translational and rota- 
tional vibrations under typical microgravity and 
terrestrial conditions for typical semiconductoi 
melts was performed. A snapshot of a typical flow 
pattern for translational vibration is presented in 
Fig. 3. Even in the total absence of gravity the 



Fie 3 Typical instantaneous 3D melt flow patterns for a lateral translational vibration at * 

R a = 0 Ra T = 7 ’S x 10 4 , Pr = 0.01. to = 100 Hz. The velocity components are V x , V y , V : , P is the pressure. T is the t mp . 

is the instantaneous temperature disturbance and \V\ is the velocity magnitude. The grayscale range corresponds to maximum va 
(white) of the velocity, temperature and pressure variables to their minimum values (black). Vibrations are applied along the tor 


(x-direction). 
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vibrations have resulted in detectable flows. For the 
cases examined, the temperature distribution le- 
mains almost unperturbed (due to the low Pr and 
weak flow strength). 

The ansle between the direction of vibration and 
the ampoule has been studied for translational vi- 
brations in the presence of an axial temperature 
gradient. At hitzh frequencies and when the angle is 
zero, no influence of the vibration on the flow w'as 
observed, even when vibrational the Rayleigh num- 
ber is verv high. The maximum observed effect 
corresponds to an angle of 90 . Here transpoit is 
significantly enhanced. 

Typical flow patterns for rotational vibrations 
flow regimes are presented in Fig. 4. Maximal 
values of flow' velocity are obseiNed at the end of 
the ampoule that is farthest from the lotation ori- 
gin. 

The influence of vibrations on heat and mass 
transfer becomes significant foi oxide melts due to 
their low' thermal diffusivity (Pr ~ 10). These flow' 
patterns are shown in Fig. 5 for the case of circular 
polarized vibration. Initially (at time t = 0), the 
species concentration was c — 1 at the lowei quai- 
ter of the cylinder and c = 0 elsewhere. The evolu- 
tion of the species concentration (process of mixing) 
is shown in Fig. 6 together with minimum and 
maximum values of velocity (for the whole domain) 
components. Complete mixing occurs in about 
10 s. The heat transfer (local Nusselt number at the 
top and the bottom) is also enhanced by about an 
order of magnitude. If the frequency of vibration is 
high, of the order of 100 Hz (for fixed Ra fi ). then the 
chances in heat and mass transfer due to vibrations 
become less significant. This corresponds to earlier 
experimental obseiwations [7.8]. 

Our results show' that both translational, circular 
polarized and angular vibration can cause average 
melt flow for a range of parameters typical of prac- 
tical semiconductor grow'th. For a given vibration 
amplitude and frequency, circular polarized and 
rotational (angular) vibrations result in moie inten- 
sive melt flows than translational ones. 

The influence of forced vibration on y-jitter-in- 
duced flow's using SAMS microacceleration data 
from the USML-2 mission (Fig. 7a) was also in- 
vestigated. Motivated by the predictions of the 
averaged equation theory presented in Ref. [11], 



Fie. 6. Temporal evolution of (a) species concentration C. min- 
imum and maximum values of C. (b) velocity extrema V v , and (c) 
V . for circularly polarized vibration applied to an oxide melt. 
Ra = 7.25 x 10 3 . Raf = 7.4 x 10 6 , Pr = 15. u> = 10 Hz. 


translational vibration was applied parallel to the 
ampoule axis (and thus, the temperature gradient) 
in an attempt to damp unwanted irregular time- 
dependent flow caused by g-jitter. While the flow 
variation with time becomes more regular, we die 
not succeed in completely suppressing the z/-jittei 
flow (see Fig. 7b). 

We found that the use of the same amplitudi 
vibration in the direction orthogonal to the am 
poule axis is more effective. This induces intensiv 
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G-jitter 



Time 



Fia. 7. Minimum and maximum values of the V x velocity com- 
ponent versus time corresponding to a melt flow response to 
g-jitter (a), g-jitter and longitudinal vibration. Ra(f = .-.2? x 10 J . 
and (b) g-jitter and lateral vibration. Raf = 7.40 x 10b 


thermal vibration flows and flow disturbances 
due to g-jitter become practically insignificant 
(Fig. 7c). 


4. Conclusions 

The influence of translational, circularly polari- 
zed and rotational (angular) vibration in analysis in 
a model Bridgman melt growth configuration was 
investigated. The nature of the flows produced by 


the types of vibration under consideration necessi- 
tated the use of the full 3D equations governing the 
transport of heat, mass and momentum. The gov- 
erning equations were solved numerically. Flow 
patterns for translational, circular polarized and 
rotational (angular) vibrations and g-jitter micro- 
accelerations were analyzed. For translational 
vibration, thermovibrational flow is strongly de- 
pendent on the angle between the vibration direc- 
tion and the temperature gradient. Circular 
polarized and rotational vibrations result in more 
intensive melt flows than translational ones. The 
simultaneous action of g-jitter and translational 
vibrations is currently being studied from the view- 
point of using applied vibration as a means of flow 
control. 
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Abstract 

A numerical solution tor thermal convection flows in a 
semiconductor melt with strong static magnetic field is 
presented. Rectangular cavity with different aspect ra- 
tios and gravity direction aligned and misaligned with 
the magnetic field vector are considered. Three numeri- 
cal methods are compared. Although, the generated flows 
have extremely low velocity, the numerical solution of the 
soverning equations involved is very complicated due to 
the thin boundary layers. It is shown that the finite ele- 
ment approach with regularization can provide numerical 
solutions in a wide range of Ha numbers (up to 10 ). The 
results compare favorably with the asymptotic theoretical 
solutions. 

1 Introduction 

Application of magnetic fields is a promising approach lor 
the reduction of convection during directional solidifica- 
tion of electrically conductive melts. Current technology 
allows the experiments with very strong static fie ids (up 
to 80 KGauss) for which nearly convection free segrega- 
tion is expected in melts exposed to stabilizing tempera- 
ture gradients (vertical Bridgman melts with bottom seed- 
ing) [1]. 

However, the reported experimental studies have 
yielded controversial results [2].[3]. Therefore, the com- 
putational methods are. a fundamental tool that may en- 
hance our understanding of the phenomena occurring dur- 
ing the solidification of semiconductor melts. Moreover, 
the effects like the bending of the isomagnetic lines, dif- 


ferent aspect ratios and misalignments between the grav - 
ity and magnetic field vectors which are difficult to model 
analytically, be studied through numerical simulations. 

The reported numerical models and results are not able 
to explain the published experimental data [4], [5]. The 
computational task is complicated because of the thin 
boundary layers [6]. although velocity of the generated 
flows is extremely low. 

Here, for comparison, three different numerical ap- 
proaches we have been used : 

(1) The spectral method implemented in [7]. 

(2) The finite element method w'ith regularization for 
boundary layers [8], 

(3) The multiquadric method, a novel method with global 
radial basis functions [9]. 

The results obtained by these three methods are pre- 
sented for a wide range of Hartmann numbers correspond- 
ing to magnetic fields B from 0.05 to 5.0 Tesla (0.5 to 
50.0 KGauss). Comparison and discussion of accuracy, 
efficiency, reliability and agreement with the asymptotic 
solution are presented. 

2 Governing equations 

The two dimensional steady state thermal convection ot 
an incompressible viscous fluid (melt) in a rectangular 
cavitv having width D and hight H was considered. More 
general geometries are discussed in [17], [18], [19]. 

The governing momentum, continuity and energy equa- 
tions (using the Boussinesq approximation) are respec- 
tive! > : 



(WIV — /VT : Y -r V/> = RaPr ■ 0 • e. - F (1) 

V-V=0 (2) 

(VT)0 = V : 0. (3) 

where length, time and velocity are scaled by L. L~ j K and 
k/L respectively. Here L is equal to the the smaller of 
the two characteristic dimensions involved ( D and H ). and 
K = i is the thermal ditt'usivity of the melt. F is a body- 
force due to magnetic field (Loren/ torce). The Lorenz 
force is given by 

F = PrHa : [{\ x e B ) x e B ) = {PrHa'V , . 0) (4) 

where V] is the horizontal component ot the velocity. The 
electrical potential was neglected, because, when a ver- 
tical magnetic field is applied, the electrical potential is 
uniform iVo — 0) [I]. The electrical potential can not 
be neelected if the symmetry is broken, e.g. it the mag- 
netic held direction and the gravity vector are slightly mis- 
aligned. However, to simplify the study, we neglect this 
effect. According to [ I ). Joule heating due to electromag- 
netic held can be neglected as well. The nondimensional 
temperature 0 is scaled by 0 — (7" — T ( , ,,,/)/ AT. A T — GL. 

G is a vertical temperature gradient r ^-. The Prandtl. 
Ravieith and Hartmann numbers are. respectively, given 
by 

v pA TgL' 

Pr = -. Ra = - • Ha = LB 0 

w vk 

Here p. e.v.p.a are thermal expansion coefficient, grav- 
itational acceleration, kinematic viscosity, density, and 
electrical conductivity: Bn is the magnetic held intensity, 
e . t b are the unit v ectors in the direction ot gravitational 
acceleration and magnetic field. 

The boundary conditions for eq. tl)-(3) are: (i) the 
non-slip condition for the velocity. V = 0. on the walls, 
and lii i specified temperature distribution. © - y on side 
walls, and © s = -0o( v-L) : at the bottom. To exclude 
the undefined constant, the pressure w as set to zero at an 
arbitrary location. /’(.Yo.yo) = 0. 

3 Problem parameters 

The problem was solved using the properties ot Germa- 
nium i Get melt exposed to magnetic fields having inten- 
sity Bn in a range from 0 to 5 Tesla land in lew cases up to 


50 Tesla). The direction of magnetic field v ector is axial. 
eg — (0. 1). The corresponding Hartmann number varies 
from Ha = 0 to 2 1 70 (and in few cases to Ha = 2. 17 • 10 ). 

The flow domain is rectangular cavity having hight H 
and width L. Two cases were considered (a) L = 1 cm. 

H = 2 cm. and (b) L = 2 cm, H = 1 cm. 

The temperature gradient on a side wall is G = 
lOK/cm. The average lateral temperature gradient at 
the bottom boundary was 0.5 K /cm. The corresponding 
Rayleigh number is Ra = 1 .24 • 10 s . Earth gravity g was 
considered to be (i) aligned to the magnetic field direction. 
g = go(0. — 1) and (ii) misaligned by 0.3 degrees. 

4 Numerical methods 

Three numerical methods were used to solve the system ol 
eq. (1) - (3) for the boundary conditions and parameters 
given above. 

4.1 Spectral element method 

To use the spectral element (SE) code NEKTON [7]. the 
domain was divided into quadrilateral elements, refined 
near the walls (Fig. la). 8 x 8 Chebyshev polynomials 
were used inside each of elements for field variables ap- 
proximation. In our tests total number of elements was 
362. The total number of unknowns was 5 10 

4.2 Finite element method with regulariza- 
tion for the Navier-Stokes equations 

The finite element method with regularization tor the 
Navier-Stokes equations (FEMR) was proposed in [8] lor 
high Re number flows. It was shown that such a regular- 
ization works well in case of flows with thin boundary lay- 
ers. even with few mesh nodes placed inside the boundary 
layer. For the present problem, the continuity equation (2) 
is modified as follows 

V-V = tV-(Vp — F — RaPr ■ © • e g ) (5) 

where x is a small regularization parameter. For t 0 eq. 
(5) approaches the original equation (2). The boundary 
condition for the pressure on the wall is 

(Vp_F-«aPr-0-tg-n = O. (b) 

where n is a unit wall normal vector. Eqs. (5) and (6) 
present the main feature of this method, and ensure the 
balance of the component of the force normal to the region 
boundary. This approach allows the use ot the same order 
finite element approximation for the velocity and pressure 




with all unknowns located at the same nodal points. The 
justification of this regularization is given in 1 101. Similar 
terms have been obtained as a result ot the consistent treat- 
ment of time-advancement for the divergence-equation by 
Lf'jhner (see [11]). Lohner has also shown that similar 
terms appear in the discrete equations as a result of diffci- 
ent order finite element approximations used tor interpo- 
lation of velocity and pressure. 

The numerical solution is insensitive to the value ot t. 

In out calculations vve have chosen the value ot i within 
the ranee 1CT 8 to 10 -4 and obtained nearly identical so- 
lutions? For a smaller value of x the discrete equations 
become nearly incompatible, and numerical solution ex- 
hibits strong spatial oscillations. 

Simple linear finite elements were used tor numerical 
approximation ot velocity, pressure and temperature. Tri- 
angular meshes with 40 x 100 and 80 x 100 nodes were 
refined near the walls (Fig. lb). Since on both meshes 
yielded nearly identical results, we used the 40 x 100 
mesh in most runs. Total number of nodes and unknowns 
was respectively 4- 1(T and 16* 10\ The FEMINA/3D 
CFD code (Finite Element Method IN Applications) [12] 
was modified to implement the proposed regularization 
method. Discrete finite element equations corresponding 
to the eqs. ( 1). (5). (3) were solved together simultane- 
ouslv by the CNSPACK solver [121 using the CGS-type 
iterative technique and high order preconditioning by in- 
complete decomposition. 

4.3 Multiquadric radial basis function 
method 

The Multiquadric Radial Basis Function (MQ) Method 
is a novel meshless collocation method with global ba- 
sis functions. The concept ot solving partial differential 
equations < PDE i using radial basis functions (RBFs) was 
introduced bv Kansa in 1990 [9]. He implemented this 
approach for the solution of hyperbolic, parabolic, and 
elliptic PDE^ using the MQ RBFs proposed by Hardy 
[ 13 ]. [ 14] for interpolation of scattered data. 

An RBF is a function that depends only upon the dis- 
tance between a point (.v.v) and a reference node (Ay.\y). 
Among studied RBFs. only the MQ RBFs have been 
proven to have an exponential convergence for the func- 
tio n interpolation [16]. A M Q RBF is given by g,{.\\y) = 
/pv _.v. 1 : - (V - . where c, is called the shape 

parameter. The numerical experiments for parabolic and 
elliptic PDEs by Kansa [9] show high accuracy and effi- 
ciency of the MQ scheme. A brief review on MQ RBF for 
the solution of PDE can be found in [15] and on the RBF- 
PDE Web site [22]. This approach results in modest size 


systems of nonlinear algebraic equations which can be ef- 
ficiently solved by using widely available library routines 
and linear solvers for dense matrices. 

For a given set of N nodes the solution tor unknow n \ . 
p or © is approximated as a sum ot MQ functions w ith the 
coefficients as unknown. These coefficients are found by 
collocating governing equations at the internal nodes and 
boundary conditions at the boundary nodes. The nonlin- 
ear algebraic system is solved by Newton method. 

.25 x 25 uniformlv distributed nodes and constant shape 
parameters cj = c*o = const were used tor all functions. 
Total number of unknowns was 2500. 

5 Results and discussion 

5.1 Convection in rectangular cavity with 
H/D=2 

The nondimensional parameters Pr = 0.006, Ra = 1 .2? ■ 

10 s we re used in all calculations. 

Flow without magnetic field: B=0, Ha=0 

The flow domain has D = 1 , H = 2 and the length scale 
was L = D. The temperature distribution at the bottom 
boundary is given by Ob = -3.575 10 (l — 4 a“). The 
results for the case a = 0 are shown in Fig. 2 (a is the 
ansle between the gravity vector and the vertical axis). 
The solution obtained by all three methods are nearlv 
identical. The flow is driven by radial temperature gra- 
dient caused by parabolic temperature profile imposed on 
the bottom boundary. The flow pattern consists of two 
counter-rotating symmetric cells, located at the lower cor- 
ners. Note that the stabilized axial temperature gradient is 
suppressing the flow. 

If the direction of gravity vector is misaligned with the 
ampoule axis by 0.5 deg., the flow pattern becomes quite 
different. The component of gravity normal to the temper- 
ature gradient becomes a driving force for the convection. 
A single roll is formed, while the magnitude of melt ve- 
locity is higher by a factor of two to three. 

Flows under magnetic fields 
5 . 1.1 B = 0.05 Tesla , Ha=21.7. 

The MQ method did not yield a solution, because the 
Newton method did not converge (since the Jacobian be- 
comes ill-conditioned). 

The solution by the SEM and FEMR methods show no- 
table differences. The SE solution for the velocity field 
exhibits numerical oscillations between the mesh nodes. 
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The how pattern from FEMR is the same as in the absence 
of a magnetic field (Fig. 3). The vertical velocity profile 
at v = 0.25 shows a boundary layer. The flow velocity is 
decreased by about a factor of two. 

5 . 1.2 B - 0.5 Tesla . Ha=217. 

The boundary layer becomes very thin, and the flow ve- 
locity is about two order of magnitude lower compared 
to B = 0 (not shown). The velocity profiles from the SE 
computation exhibit spatial oscillation with velocity sign 
chanae between mesh nodes. The FEMR can provide the 
results still without difficulty, the velocity profiles remain 
smooth. 

5 . 1.3 B =5.0 Tesla, Ha=2 170. 

The results from the SE computation showed strong nu- 
merical instability. The FEMR solution indicates that the 
flow' pattern is about the same, while the velocity i^> lower 
about two order of magnitude compared to B = 0.5 Tesla. 
The boundary layers are extremely thin (0.0 lev?:*. and 
therefore almost invisible on a plot (Fig.4). 

In case of a misalignment of gravitational acceleration 
with ampoule axis, the flow pattern changes to one big cell 
for this and all other values of magnetic field considered. 

5.1,4 B = 50 Tesla , Ha=21700. 

This case was done just to test the ability of proposed 
FEMR method, the solution still remains smootr. with 
even three times more thin boundary layer compared to 
B = 5.0 Tesla. Stream lines for this case are shown in 
Fig. 5. and the velocity amplitude is presented by the most 
right point on a plot in Fig. 6. 

Stretching of the stream lines caused by the magnetic 
field is shown in Fig. 5. This stretching pre>ented 
schematicalh in [1]. [5] but. to our knowledge, compu- 
tational results were not shown. 

Discussion 

Figure b shows the maximum radial velocity calculated, 
using the FEMR method, for different values of me im- 
posed magnetic field B. The maximum ot horizontal (ra- 
dial) velocity \ersus B is presented by few curves, marked 
as "Vr(bY* for FEMR on 40x100 mesh uniform in the 
vertical direction, by "V 7 r{f) for 40x100 mesh renned at 
the walls and bv "Vr{d) tor 80x100 mesh refine^ at the 
walls. Results for misaligned case are presented by the 
curve w ith squares, labeled as "Vr(a = 0.5) . One can ob- 
serve the predicted asymptotic dependence V max - Ha 


[11 for all the cases, starting at about B — 0.05 Tesla 
(Ha « 20). 

The main computational difficulty of this problem is 
due to the viscous flow with thin boundary layer. De- 
spite the fact that actual flow' velocities are very low and 
the Reynolds number obtained using the computed ve- 
locities. is Re ~ 10' 1 to 10~ 6 , a big value of the Hart- 
mann number results in a relatively small coefficient at the 
highest derivative of the velocity in the momentum equa- 
tion., Solution of such a problem exhibits thin boundary 
layer with the thickness 5 ~ Ha~ \ and the the equiv- 
alent” Reynolds number Re eqv ~ Har , for B=0.5 Tesla 
Re eqv = 4.1 ■ 10 4 and B=50.0 Tesla Re eqv = 4.7 - 10 . 

5,2 Thermal convection in rectangular cav- 
ity with aspect ratio H/D=0.5 

The nondimensional temperature distribution at the bot- 
tom boundary is given by © b = -7.150- 10 (1 - a~)- 

The axial temperature gradient applied on the vertical wall 
is a i S o G = 70 KVcm. 

Flow without magnetic field: B=0, Ha=0 

The solution obtained by all three methods are close to 
each other. The flow pattern consists of two counter- 
rotating symmetric cells, that occupy most of the volume. 
Fig. 7. 

In the case of the gravity misalignment with the am- 
poule axis direction by 0.5 deg .. the axial temperature gra- 
dient becomes a main driving force for the thermal con- 
vection. This results in the change of flow pattern that 
becomes consisting of one big convective cell. 

Flow' under magnetic fields 

The results are shown in Fig. 8. Again when Ha number 
is high, all the methods except FEMR, exhibit the same 
difficulties as in a case of aspect ratio H/D = 2. A sum- 
mary of the results is shown in Fig. 9. The suppression ot 
the flow is essentially same efficient as before with sim- 
ilar asymptotic dependences — Ha~~. The velocitx 
profile in the boundary layer obtained by FEMR is shown 
in Fig. 10. One of the advantages of FEMR is that its so- 
lution remains smooth even at the big change ot the slope. 
One can see that the thickness of the vertical boundary 
laver is in agreement with asymptotic solution, 5 ~ Ha 
The tangent velocity derivative at the boundary decreases 

with Ha number as ~ ^ ~ Ha " 1 . 

Comparing between Fig. 6 and 9 it is found that mis- 
alignment's impact on the reducing of the convection is 
more important for the aspect ratio H /D - 1 . 
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Conclusions 

We compared three different numerical methods tor the 
solution of thermal convection flows in a semiconduc- 
tor melt with strong static magnetic held applied. These 
are spectral element method, finite element method with 
regularization for the Navier-Siok.es equations and mul- 
tiquadric method, a method with global basis [unctions. 
Although the generated flows are extremely low. the com- 
putational task is very complicated because of the thin 
boundarv layer at high Hartmann numbers, Ha ^ 1. We 
considered melt region geometry with different aspect ra- 
tios. and gravity direction aligned and misaligned with the 
magnetic held vector. The comparison shows that the fi- 
nite element approach with regularization can obtain sta- 
ble and reliable solutions in a wide range ot Ha number, 
up to 10 4 . These results compare favorably with asymp- 
totic solutions. 

The main difficulty of this problem is that a flow has 
a verv thin boundary layer. Despite the tact that ac- 
tual Revnolds number is very low. Re ~ 10“ 1 to 10 6 

a hish value of the Hartmann number results in a rel- 
ativelv small coefficient at the velocity Laplacian in the 
momentum equation. Solution ot such problem exhibit 
thin boundary layers with related, like tor high Reynolds 
number flows, difficulties. That is one of the reasons tor 
the discrepancy in the results that numerical studies re- 
ported. Both the spectral method and the multiquadric 
method with global basis functions needs improvement to 
deal with thin boundary layers. Multilevel approximation 
bv Fasshauer [20].[21] can be one of the ways. 

Numerical solution ot these problems by available 
commercial CFD codes may be not efficient or not possi- 
ble. Adaptive algorithms can be a promising solution. De- 
velopment of more accurate and efficient solution meth- 
ods for this problem is necessary. 
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Fiaure 3: Thermal convection with magnetic held tor geometry 1. D = 1. H = 2. Ra — 1.25 • 1 . a - _ . ( 
0.05 Tesla): tat stream function for aligned and ib) misaligned by 0.5 degree gravity direction, (c) vertical ve oci \ 

profile V,.(.v)/(0.225 cm/s) for (a) at y/H = 0.25 by FEMR. 



Figure 4: Thermal convection with magnetic field for D - 1 . H - 2. Ra - 1 .25 • 10- . Ha 2170. 0 (B 5 d) 7>5/a I 

(a ; stream function: ib) vertical xelocity profile V, v)/(0.225 cm/s. at y = 0.25: tc) velocity profile 1 >a jc)/( 0.^5 cm/s 
tor misaligned b\ 0.5 deszree izravity direction at > 






Figure 5 : Strcatching of the (low streamlines with increasing magnetic field, (a) B 
, ^correspond. respectively, to B = 0.05.0.5.5.0 and 50.0 Tesla 


0. <b) B 


0.005 Tesla, and (c) to 
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VELOCITY VS MAGFIElD H=2.D*f 



Figure 6: A summary of magnetic field suppresion of the flow for H = 2. D = 1. Maximum value of horizontal (radial) 
velocity versus B: Vr(b) on 40x100 uniform in vertical direction mesh. Vr(/) tor 40x 100 refined near walls mesh, 
and Vr(d). 80x100 mesh refined at the walls. Predicted asymptotic dependence V„ iav ~ Ha - is observed tor all the 
cases including the misaligned one. starting from about B = 0.25 Tesla. 



Future Thermal convection without magnetic field for D = 2. H = 1 . Ra = 1.25- 10" : (a) stream function, (b) stream 
function for misaligned configuration, tc) vertical velocity profiles by SEM and FEMR. V, (.v) tor (a) at v - 0-5. and 
,di vertical velocity profiles by SEM and FEMR. V v (.v) for (b) at y = 0.25. Velocity scale is 0.225 cm/s. 
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Figure 8: Thermal convection with magnetic field. D = 2. H = \. Ra = 1.25- 10\ Ha = 2170 ( B = 5.0 Tesla): (a) 
stream function, (b) stream function for misaligned case. (c) vertical velocity profile \\ (a) tor (a) at y/H = 0.25 and 
(d) vertical velocity profile Vyix) tor tb) at y/H — 0.25: velocity scale is 0.225 cm/s. 



Figure 9: Summary of magnetic field suppresion 
dial) velocity versus B for aligned and misaligned 
observed for both cases. 


of the flow for H = 1. D = 2: maximum value of horizontal (re- 
configurations. Predicted asymptotic dependence V mcLX — Ha~- is 
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Introduction 

The character of natural buoyant convection in rigidly contained inhomogeneous fluids can be 
drastically altered by vibration of the container. For certain experiments and operating 
conditions, vibrations are expected to have a significant influence on heat and mass trans er 
onboard the International Space Station (ISS). Furthermore, it appears that g-jitter vibrations will 
exist on ISS over a wide range of frequencies [1]. In general, vibrational flows are very comp ex 
and are governed by many parameters. This complexity makes it almost impossible to correctly 
predict vibrational effects empirically. Thus, a careful theoretical approach combine wit 
numerical modeling is essential. Available flight experiment data clearly show that, once 
initiated by “g-jitter", the effects of convective flows can persist for long times even when the g- 
jitter disturbance (and consequent flow') were short-lived [2-7]. 

In many terrestrial crystal growth situations, convective transport of heat and constituent 
components is dominated by buoyancy driven convection. Control of convective transport 
continues to be an important aspect of crystal growth research. Several groups are y 

pursuing control of convection using static and rotating magnetic fields. Magnetic fields cannot 
be used for flow control in melts and solutions that are poor conductors. Flow suppression 
through vibration or vibro-convective mixing may offer an attractive alternative in such cases. 

Recent work has shown that the character of natural buoyant convection in non-uniformly heated, 
rigidly contained inhomogeneous fluids can be drastically altered by vibration of the container. / 
review and relevant theoretical and experimental research can be found in publications [ - ]• 

Thus. vibrational induced flow can potentially be used to influence and even control transport in 
some crystal growth situations. A practical quantitative understanding of vibrational convection 
as a control parameter in crystal growth situations is currently not available. The objective of t e 
work is to assess the feasibility of the use of vibration to suppress, or control, convection in order 
to achieve transport control during crystal growth. 

Buovancy driven vibro-convective motion occurs when oscillatory displacement of a containei 
wall induces the acceleration of a container wall relative to the inner fluid. The vibration mav e 
viewed as a time-dependent modulation of steady gravity. In a closed container the fluid wi 
move as a rigid body with a container. If. however, the fluid density is nonuniform, fluid motion 
mav ensue. The magnitude of this motion, of course, depends on the orientation of the vibrationa 
direction with respect to the local density gradients. Note that, similar to Rayleigh-Benard 
configurations, there mav be a -critical" threshold for the coupled vibrational frequency and 
amplitude, to cause connection. Interestingly, it should be noted that in case of a constant density 
fluid subject to spatially nonuniform vibration, fluid motion can also occur (for example, angular 

vibration [11]). 

To properly investigate influence of translational, circularly polarized and rotational (angular) 
vibration necessitates the use of the full 3D equations governing the transport of heat, mass and 
momentum. Selected examples of our ongoing work on this topic are outlined below . 


We consider a purely thermo-vibrational convection in a differentially heated cylindrical cavity 
with no consideration of solidification. The fluid is taken to be Newtonian and the Boussinesq 
approximation is assumed to hold. The calculations were performed tor identification and 
characterization of thermovibrational flow and are part of an ongoing project invoh in = t o\\ 
visualization model experiments being conducted by Feigelson [10]. 

Mathematical model for translational and polarized vibration 

Translational vibration corresponds to a linear displacement such as, for example, u-dcoscot. 
where d is a real vector giving the displacement magnitude and co is the frequency. In this case 
the ampoule is displaced back and forth upon the same line. Polarized vibrations are 
characterized by a displacement u = Re{de M \ where d = d y - id 2 (see Fig. l(a.b)). Here the 

instantaneous vibration direction rotates in the polarization plane defined by the real vectors 
d i and d ~ . A sketch showing both translational and circular polarized vibrations is presented in 
Fis. 1 . In a reference frame fixed to a vibrating ampoule, the momentum equation is 


+ (v-V)v = -Vp + PrV : v + Ra T Pr{0 + aC)k + Ra T Pr{0 + aC)/(£2.r), (1) 


where length, time and velocity are scaled by Rq, Ro 2 /Kand k/R 0 . Here Ro is the ampoule radius 
and K is the thermal diffusivitv. The nondimensional concentration and temperature, are given by 
0 . and C. respectively. The function f{Q.t) is the acceleration of the vibrating ampoule and 

Q =cu/? ( ; / k is a dimensionless frequency. The continuity and heat-mass transfer equations 
complete the problem formulation. The Prandtl. Schmidt, thermal and vibrational Rayleigh 
numbers and the buoyancy ratio are given by 


Pr = —, Sc = 

K 


V 

D 


Rcij - 


MML a= M^,Ra* PclAHl 


VK 


(5 AT 


VK 


( 2 ) 


Here and j3 c are the thermal and solutal expansion coefficients and AT. c m , g, d. co. k. v, D are 
the characteristic longitudinal temperature difference, reference concentration in the melt 
gravitational acceleration, vibrational displacement amplitude and frequency, direction of 
gravity, kinematic viscosity and solute diffusivity, respectively. The dimensionless number Ra* T 
is the vibrational Rayleigh number and Ra* s = a Ra* T . Equation (2) is solved together with the 
equations governing heat and species transfer and the continuity equation. 

Rotational vibration 

The equations of motion for angular vibrations take on a more complicated torm (see Fig. lc). A 
container of length L is subjected to an angular displacement 0(r) in the plane. Here the 

coordinates .v -are referred to a fixed laboratory' frame of reference. The position vector to the^ 
mass center of the cylinder is parallel to the side of the cylinder and is given by q - Ro[-sin6\\ + 
cos $'.*] where R 0 is the distance from the origin 0* to the mass center of the cylinder and 6{t) = 
esinQnt. In a frame of reference moving with the container, the equations of motion have the 
form 


( 3 ) 


p£L = JjvT - pk ( sin 0/, + mv0i\! + p[2i2v - Q'X + &x + £'*2 0 ‘^o cosQ a ti^ 

Dr L 

vvhere D/Dt denotes the material derivative, r is the velocity of the fluid relative to the movin 
reference frame, p is the density of the fluid. Q is the rate of rotation tensor for the moving 

frame with respect to the fixed frame of reference. Q is its time derivative and T is the 
Newtonian stress tensor for the fluid. The dimensionless equations governing the transport ot 
momentum, mass and heat in the cylinder are obtained after using L, L 2 /k, k/L, and A h - c 
to scale respectively, length, time, velocity and temperature. The governing dimensionless 
parameters are the dimensionless frequency Q=Q G L 2 /k, , the dimensionless container radius, 0 
= Rq/L. the Prandtl. Pr. and the Rayleigh. Ra. vibrational Rayleigh Roq. and Ekman. E. 
numbers. The latter are given by 



pSTgL_ 
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, Ra a 


£= v 

kv El,, Is 


(4) 


where B v q and K are the coefficient of thermal expansion, kinematic viscosity, gravitational 
acceleration and thermal diffusivity. respectively. This system of equations differ from the usual 
equations in the absence of rotation in that additional terms are present: the Coriolis term which 
is proportional to ePr/E. and the centrifugal term which is proportional to £ Wa<pPr and varies 
with linearlv with position in the ampoule. The importance of the latter term depends on the 
dimensions of the amplitude of the angular vibration, £. and the ratio 0. The rocking motion ol 
the angular vibration under consideration means that centrifugal terms give rise to a periodic 
forcing that fluctuates about the mean value at twice the period ot the angular vibration. 


Since the above svstem of equations has not been well studied, a conservative approach was 
adopted for the studv of angular vibrations and we confine our investigation to a P ar amemc stu > 
of flow regimes and'transitions for thermo-vibrational situations in the absence of solidification. 


Solution method 

The equations are solved in primitive variable form using a Finite Element Method code 
FEMINACD developed bv the PI [14]. The continuity equation and momentum equations are 
considered simultaneously at each time step. This eliminates many problems related to boundary 
conditions and places only slight limitations on the time step size for transient problems^ The 
regularization for the incompressibility condition makes the solution procedure more e 1 - 

and allows the same order finite element approximation for both the velocity and P r ^sure [15], 
This approach makes it possible to solve large 3D time-dependent problems (up to ,00.000 
unknowns! on a SGI 02 workstation with reasonable computation times. 

We implemented the above 3D models of convective buoyancy-driven melt flow in differentially 
heated cylindrical containers using the FEMINA/3D code. This code was carefully tested on 
benchmarked experimental and numerical data for a variety of 2D/3D viscous an t ermo 
convective flow problems and flows under magnetic field [15.16,19]. 

For rotational vibrations the Ekman number can be of the order 10' 4 to 10‘ ? for frequencies on^ 
the order of 1Hz. This results in large coefficients. Pr/E. for the Coriolis terms in the go\erm - 



equations anti causes difficulties in the numerical solution. To resolve this we implemented a 
hh!h accuracy solution method using preconditioning by high order incomplete decomposition 

( accuracy up to 1 O'"). The computation times reduced by one to two orders of magnitude and the 
memory size hr a factor of 8 for 3D flows compared to currently available commercial codes 
(e.2. CFD2000). A typical solution time for a transient problem is about two hours on a SO - 

Results and discussion 

We verified the validity of the Boussinesq model for semiconductor and oxide melts under 
microsravitv conditions. This topic has been discussed recently by Perera. Sekerka [17] 
Pukhnachev [18] and Gershuni. Lyubimov [11]. If the nondimensional criteria, propose y 
Pukhnachev Pit = glVV' is less than 1. then the Boussinesq model for thermal convection 
mav not be valid. Our estimates show that the Boussinesq model is quite adequate for a 
differentially heated closed ampoule and the range of parameters and material properties un er 
investigation. The values of Pit are of the order 10 4 to l(k for semiconductor and oxide melts for 
2/ go ="l0' ? to lCf\ clearly well above 1 . 

A parametric studv of translational and rotational vibrations under typical microgravity and 
teireslrial conditions for typical semiconductor melts was performed. A snapshot of a typical 
How pattern for translational vibration is presented in Fig. 2(a). Even m the total absence ot 
gravity the vibrations have resulted in detectable flows. For the cases examined, the temperature 
distribution remains almost unperturbed I due to the low Pr and weak flow strength). 

The angle between the direction of vibration and the ampoule has been studted for translational 
vibrations in the presence of an axial temperature gradient. At high frequencies and when e 
rS” zero, no influence of the vibration on the flow was observed, even when vtbta.tonal the 
Raleigh number is very high. The maximum observed effect corresponds to an angle of 90 
de°rees Here transport is significantly enhanced. 

Typical flow patterns for rotational vibrations flow regimes are presented in Fig. 2(b). Maximal 
velocity values are observed at the end of the ampoule that is farthest from the rotation ortgtn. 
The influence of vibrations on heat and mass transfer becomes significant for oxide mielts due to 
their low thermal diffusivit v t Pi - 10). These flow patterns are shown in Fig. 3(a) tathecase o 
circular polarized vibration. Initially tat time / = 0), the species concentration was c- at the 
lower quarter of the cylinder and c = 0 elsewhere. The evolution o the speces concentrmton 
(process of mixing) and velocity (minimal and maximal values of V : ) is shown in Ft c . 3(b. ). 
cTn p ete n xing occurs in about ten seconds. The hea, transfer (local Nusselt number a. the top 
and the bottom) Is also enhanced by about an order of magnitude. If the frequency of v.brauon ,s 
hurt, Of the order of 100Hz (for fixed Run), then the changes in heat and mass transfer due to 
vibrations become less significant. This corresponds to earlier experimental observations [7.8], 

Our results show lhat both translational, circular polarized and angular vibration can cause 
average melt flow for a range of parameters typical of practical semiconductor growt i or a 
given'vibration amplitude and frequency, circular polarized and rotational (angular) vtbrattons 
result in more intensive melt flows than translational ones. 

The influence of forced vibration on g-jitter induced flows using SAMS microacceleration data 
from the USML-2 mission w as also investigated [13]. Motivated by the predtettons ot the 
aveta«ed equation theory presented in Ref. [11], translational vibration was applied parallel to 
the ampoule axis land thus, the temperature gradient) in an attempt to damp unwanted tnegu ar 


time-dependent flow caused by g-jitter. While the flow variat.on with time becomes ; more 
regular we did not succeed in completely suppressing the g-jitter flow. We found that the use of 
the same amplitude vibration in the direction orthogonal to the ampoule axis is more effec i\e. 
This induces intensive thermal vibration flows and flow disturbances due to g-jitter become 

practically insignificant. 

Summarv 

The influence of translational, circularly polarized and rotational (angular) vibration m analysis 
in a model Bridgman melt growth configuration was investigated. The nature of t e ows 
produced bv the types of vibration under consideration necessitated the use of the tu 
equations governing the transport of heat, mass and momentum. The governing equations were 
solved numerically. Flow patterns for translational, circular polarized and rotational (angu ar) 
vibrations and g-jitter microaccelerations were analyzed. For translational vibration, 
thermovibrational flow is strongly dependent on the angle between the vibration direction and 
temperature gradient. Circular polarized and rotational vibrations result in more intensive melt 
flows than translational ones. The simultaneous action of vibrations and magnetic field [19] is 
currently being studied. 
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Fig. l:Translational vibration (a), d, or d 2 = 0 : polarized vibration (b), d h d 2 * 0 ; cp is the angle between gravity 
vector and the ampoule axis, a is the angle between the vibration direction and the ampoule axis: (c) rotational 
(angular) vibration. The container is rotated at an angle 9(t) about a center of rotation at x* = 0 . The vector q* 
connects the center of rotation to the mass center of the container. 
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(b) 


Fis. 2: (a) Instantaneous 3D flow patterns for a lateral translational vibration at Og. Ra = 0. Ra T - 125 -10 . 
Pr=0.01 . co = 100 Hz. The velocity components are V„V„V.. P is the pressure. T is the temperature. AT is the 
temperature disturbance and IV1 - velocity magnitude. The grayscale range corresponds to maximum values (white) 
of the velocity, temperature and pressure variables to their minimum values (black). Vibrations are applied along the 
horizontal (x-direction); (b). 3D melt flow patterns for angular vibration at zero-g, Ra=0, Ra a - 4.6 ■ 10 , Pr -0.01. 


CO = 100 Hz. 



Fia. 3: (a) 3D flow patterns for circularly polarized vibration, Ra - 7.25 • /0\ Rci 7.4-10 , Pr 15, co 10Hz., C is 
the concentration: (b) Temporal evolution of velocity extrema V., and. (b) species concentration C. 
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The Multiquadric Radial Basis Function (MQ) Method is a meshless collocation method with 
global basis functions. It is known to have exponential convergence for interpolation prob- 
lems We descretize nonlinear elliptic PDEs by the MQ method. This results in modest-size 
systems of nonlinear algebraic equations which can be efficiently continued by standard contin- 
uation software such as AUTO and CONTENT. Examples are given of detection of bifurcations 
in ID and 2D PDEs. These examples show high accuracy with small number of unknowns, as 
compared with known results from the literature. 


1. Introduction 

Nonlinear multidimensional elliptic partial differen- 
tial equations (PDEs) are the basis for many sci- 
entific and engineering problems, such as viscous 
fluid flow phenomena, chemical reactions, crystal 
growth processes, pattern formation in biology, etc. 
In these problems it is crucial to understand the 
qualitative dependence of the solution on the prob- 
lem parameters. 

During the past two decades the numerical con- 
tinuation approach has become popular for qual- 
itative study of solutions to nonlinear equations, 


see e.g. [Rheinboldt, 1986; Allgower k Georg, 
1990; Doedel et ai , 1991; Seydel, 1998] and ref- 
erences therein. Several software packages, such 
as AUT097 [Doedel et ai , 1997] and CONTENT 
[Kuznetsov k Levitin, 1995-1997], are currently 
available for bifurcation analysis of systems of non- 
linear algebraic equations and ODEs, with only 
limited bifurcation analysis for ID elliptic PDEs. 
For 2D PDEs, we mention the software package 
PLTMG [Bank, 1998] that allows to solve a class of 
boundary value problems on regions in the plane, 
to continue the solution with respect to a pa- 
rameter and even to compute limit and branching 
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points. This software combines a sophisticated fi- 
nite element discretization with advanced linear al- 
gebra techniques. Numerical continuation for ID 
and 2D elliptic PDEs is currently an active re- 
search area, see e.g. [Neubert, 1993; Shroff &: Keller, 
1993; Schwetlick et al., 1996; Chien et al., 1997; 
Davidson, 1997; Kuznetsov et al., 1996; Chien 
& Chen, 1998; Doedel & Sharifi, to appear] and 
[Govaerts, 2000, Chap. 10] for reaction-diffusion 
equations; and [Mamun & Tuckermann, 1995, 
Poliashenko & Aidun, 1995] for CFD. The typi- 
cal approaches used are based on the finite ele- 
ment or finite difference discretization of the PDEs. 
They result in very large (thousands or tens of 
thousands for 2D problems) systems of nonlinear 
algebraic equations with sparse matrices. The 
continuation process is typically based on the 
predictor-corrector algorithms that require solving 
nonlinear systems by the Newton type method at 
each continuation step. For the bifurcation analysis 
during the continuation process, one usually needs 
to compute at least a few eigenvalues of the Jaco- 
bian matrix at each continuation step. The meth- 
ods currently used both for the continuation and 
the corresponding eigenvalue problems are variants 
of Krylov subspace methods and recursive projec- 
tion (RPM). Solving the resulting linear system and 
the eigenvalue problem require sophisticated algo- 
rithms and considerable computer resources (CPU 
time, memory, disk space, etc.). 

In this paper we report results of numerical 
experiments with continuation and detection of bi- 
furcations for ID and 2D elliptic PDEs discretized 
by the Multiquadric Radial Basis Function (MQ 
RBF or, simply, MQ) method. The MQ method, in- 
troduced for solving PDEs in [Kansa, 1990a, 1990b], 
is a meshless collocation method with global basis 
functions which leads to finite-dimensional prob- 
lems with full matrices. It was shown to give 
very high accuracy with a relatively small number 
of unknowns (tens or hundreds for 2D problems). 
The corresponding linear systems can be efficiently 
solved by direct methods. This opens a possibil- 
ity for using standard continuation software, such 
as AUTO and CONTENT, designed for bifurcation 
analysis of modest-size problems. We also note that 
the MQ method does not require predetermined lo- 
cation of the nodes as the spectral method does. 

In Sec. 2 we summarize previous results on 
solving PDEs by the MQ method. In Sec. 3 we 
formulate an adaptation of the MQ method suit- 
able for the discretization of parametrized elliptic 
PDEs. In Sec. 4 we present results of our numerical 


experiments with a ID eigenvalue problem and in 
Sec. 5 we present results of our numerical experi- 
ments with continuation of solutions and detection 
of bifurcations for ID and 2D elliptic PDEs. In 
Sec. 6 we discuss our results. 

2. Review of the MQ Method 
for Elliptic PDEs 

The concept of solving PDEs using radial basis 
functions (RBFs) was introduced by Kansa in 1990 
[1990a, 1990b]. He implemented this approach for 
the solution of hyperbolic, parabolic, and ellip- 
tic PDEs using the MQ RBFs proposed by Hardy 
[1971, 1990] for interpolation of scattered data. 

There exists an infinite class of RBFs. An 
RBF is a function /j(x) & R, x € R (say, in 
ID case) that depends only upon the distance be- 
tween x and a reference node Xj. A MQ RBF is 
9j ( Cj , x ) = ((x-Xj ) 2 +c)) x < 2 , where Cj is called the 
shape parameter. In a comprehensive study, Franke 
[1982] compared (global) RBFs against many pop- 
ular compactly supported schemes for 2D interpo- 
lation, and he found that the global RBF schemes 
were superior on six criteria. Among the studied 
RBFs still only the MQ RBFs are proven to have 
an exponential convergence for the function inter- 
polation [Madych & Nelson, 1990; Wu & Shaback. 
1993], Madych [1992] showed theoretically that 
the MQ interpolation scheme converges faster as 
the constant MQ shape parameter becomes progres- 
sively larger. 

The numerical experiments for parabolic and 
elliptic PDEs by Kansa [1990a, 1990b] and 

Golberg and Chen [1996] show high accuracy and 
efficiency of the MQ scheme. For example, for 
a ID convection-diffusion problem, Kansa [1990b] 
showed that the MQ solution with 20 nodes had the 
maximum norm error within 10 -4 , while a second- 
order finite difference scheme with K = 200 nodes 
and an optimal combination of the central and up- 
wind differences for the problem resulted in a much 
larger error of 3 • 10 -2 . 

In the numerical experiments with a nonlin- 
ear time-dependent problem modeling the ID von 
Neumann blast wave Kansa [1990b] compared the 
exact solution and its derivatives with the MQ so- 
lution (35 nodes) and with a second-order finite dif- 
ference one. The error in the maximum norm for 
pressure, density, energy and their gradients was 
10" 6 or less for the MQ method, and in the range 
from 10~ 4 to 10 -2 for the finite difference method 
with 5000 nodes. 


Golberg and Chen [1997] showed that the solu- 
tion of the 3D Poisson equation in an ellipsoid could 
be obtained with only 60 randomly distributed 
nodes to the same degree of accuracy as a FEM 
solution with 71 000 linear elements. 

Sharan et al. [1997] showed that the MQ 
method yields accurate solutions for 2D Poisson and 
biharmonic equation, and that the MQ approach 
is simple to implement on domains with irregular 
boundaries. Cook et al. [1993] noted many benefits 
of using MQ RBFs to solve an initial value problem 
for a 3D nonlinear equation for the collision of two 
black holes. The resulting discrete system had 2000 
unknowns and was solved directly. 

Buhmann [1995] showed that RBFs and, in 
particular, MQ RBFs are useful for construct- 
ing prewavelets and wavelets. Wavelets are most 
frequently used in time-series analysis, but there 
are results for using wavelets to solve PDEs 
[Fasshauer & Jerome, 1999; Narcowich et al., 
1999]. As Buhmann points out, one can gener- 
ate true wavelets by an orthonormalization process. 
The wavelets are an elegant way to achieve the 
same results as rnultigrid schemes. The MQ RBFs 
are attractive for prewavelet construction due to 
exceptional rates of convergence and their infinite 
differentiability. 

The paper by Franke and Schaback [1998] pro- 
vides the first theoretical analysis for solving PDEs 
by collocation using the RBF methods. 

Kansa and Hon [1998] studied several methods 
for solving linear systems that arise from the MQ 
collocation problems. They studied the 2D Pois- 
son equation, and showed that ill-conditioning of 
the system could be circumvented by using block- 
partitioning methods. 

Kansa [1990b] introduced the concept of vari- 
able shape parameters Cj in the MQ scheme that ap- 
peared to work well in some cases. In [Kansa & Hon, 
199S], a recipe for selecting c 3 based upon the local 
radius of curvature of the solution surface was found 
to give better results than a constant Cj MQ scheme. 
Kansa and Hon [1998] tested the MQ method for 
the 2D Poisson equation with a set of exact solu- 
tions F = exp(ax -I- by), cos (ax 4 by), sin(ax 4 by), 
log(ax 4 by 4 c), exp(-a(i - 1/2)2 - b(y - 1/2)2) 
and arctan(ax 4 by). They obtained an accuracy 
up to 10 -5 using a modest size, 121. set of nodes, 
while locally adapting the shape parameter c y 

The multizone method of Wong et al. [1999] 
is yet another alternative method for improving 
computational efficiency. This method is readily 
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parallelizable, and the conditioning of the resulting 
matrices are much better. 

Hon and Mao [1998] showed that an adaptive 
algorithm that adjusted the nodes to follow the 
peak of the shock wave can produce accurate re- 
sults in ID Burgers equation with only 10 nodes, 
even for steep shocks with Re = 10 4 . 


3: Discretization of Nonlinear Elliptic 

PDEs by the MQ Method 


We consider the second-order system of n para- 
metrized nonlinear elliptic partial differential 
equations 


D(a)Au - /(Vu, u, x, y, a) = 0 , 
«eR, u(-), /(•) er, (x, y) e 0 c R 2 , 


( 1 ) 


where D(a) is a positive diagonal nxn matrix that 
is dependent smoothly on a, subject to boundary 
conditions 




i 


du 




u, x , y, a 


3 


= 0. 


f b (-) 6 R n 


rot 


Here a is a control parameter, and we are interested 
in studying the dependence of the solutions to the 
boundary-value problem (1), (2) on a. 

We discretize the continuous problem by the 
multiquadric radial basis function (MQ) method 
[Kansa, 1990a, 1990b; Madych &z Nelson, 1990] as 
follows. Introduce a set ©h of nodes ( N internal 
and Nb on the boundary) 


Oh = {(*«> yi)|t-i,JV C O, 

(x u yi)\i=N+l,N+N b c 90} 

and look for the approximate solution to (1), (2) in 
the form [Madych & Nelson, 1990] 

j=N + N b 

xt/i (nr, y) = ao + ^ ] a j9ji c j > 9 ) ■ 

j=i 


where '£ j Zi +Nb a i = °- We use this relationshi P to 
eliminate ayv from (3) which results in 


Uh(x, V ) 


j=N - 1 

= a 0 + x ’ - 9n{cn . x , y )) 

i=i 


j=N+N b 

+ 5Z a j(9j( c j- x ' y) ~ 9 n{ c N> X, y)) 1 

j=N+\ 

(4) 
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where dj S R n are the unknown expansion coeffi- 
cients and 


this end, we first eliminate a 2 from (4) by sub 
stituting (7) into (4) to obtain 


9 ]{Cj, x, y) = yj { I - x 3 ) 2 + (y - y } ) 2 4- cj , 
j = 1 , . . . , N + N b . 

are the MQ basis functions, and Cj > 0 are called 
shape parameters [Kansa, 1990b]. We then substi- 
tute u/,(x, y) into (1), (2) and use collocation at the 
nodes &h to obtain a finite-dimensional system 


y) 

j=N - 1 

= a 0 + a ]( 9 j(cj, x, y) - 9n{cn, x, y)) 

j= i 

j=N+N b 

+ ( a ’ 1 j > x ! y ) 

j=N+l 


ip x {a\ a 2 , a) 

= D{a)Au h {x t , y { ) 

- /(Vu /l (x,, y,), Uh(x t , j/*), x,, yi, a) = 0, 

* = 1 iV, (5) 


• a ”> Q ) 


= / 6 ( 9u ^ ,y< - -, Vi). *«. yi, a) =0, 


i — N + 1 , . . . , N ■+■ N b , 


( 6 ) 


where a 1 = (ao, . . . , ajv— l) G R n/ '”, — 

(a,v+i, a^vj € R nx ' Vb . We next modify the 


discretized system to make it more suitable for con- 
tinuation and bifurcation analysis. (1) We elimi- 
nate a 2 , associated with the boundary nodes, so as 
to minimize the number of unknowns. (2) We refor- 
mulate the resulting problem in terms of (internal) 

nodal values u, = u^(x,, yi), t = 1 N, so as to 

have the correct eigenvalue problem (to avoid deal- 
ing with matrix pencils) for the Jacobian matrix of 
(5) for detecting bifurcations during the continua- 
tion process. This is accomplished as follows. 


1. We solve (6) for a 2 (assuming that the Implicit 
Function Theorem is applicable here) to obtain 


a 2 = ^(a 1 , a), or, in components, 

dj = ipj(a}, a) , (7) 

j — N + 1 , . . . , N + N b . 


Substituting this into (5) and using the notation 
<~P — (y i , ■ • • , ¥hv), yields 

^(a 1 . [/-’(a 1 , a), a) = 0. yd-) € K' IX ' V . (8) 


2. We now want to reformulate (8) in terms of the 
nodal values U = (uj, Ut... . . u_\) € R nXiV . To 


- 9n( c n > x, y)) . 

(9) 

We now define the map T : a 1 i — > U = T(a 1 ). 
For i — 1, . . . , N: 

j=N - 1 

U{ = a o + ^ ( 9]{ c j > x i< y») 

j=i 

— 9n{ c n , Xj, yi))Qj 
j=N+N h 

+ 5Z Xi, Vi) 

j=N + 1 

- yjv(c,v, Xi, y t ))^j(a 1 , a) (10) 

Finally, substituting a 1 = T _1 (f/) into (8), 

we arrive at the finite-dimensional continuation 
problem 

G{U, a) = 0 , U, G{-) € R nx/V , a e R, (11) 


where 

G(C7, a) = tp{T a), a) , 
r : -> R jV , ip{-) € R Tlx,Vb • 


Remark 1. Note that in the case that the bound- 
ary condition (2) is linear, Tpj are linear , and con- 
sequently T is an N x N matrix. 

In Sec. 5 we consider examples of continuation 
of ID PDEs with Q = (0, 1) and 2D PDEs with 
n = (0, 1) x (0, 1). In all 2D examples we have the 
same number of nodes N s in x and y directions. We 
choose a constant shape parameter c 3 = s/{N s - 1). 
Our typical choice for s is 4 < s < 12. 

We use two types of node distributions. In 
the case of uniform node distribution (xfc, yi) = 
(kh. lh ), k, l = 0 ,..., N s , h = 1 /N,. In the case 



of nonuniform node distribution, the nodes adja- 
cent to the boundary are placed at the distance 
h = h x h from dtt, 0.1 < h x < 0.5, while the remain- 
ing nodes are distributed uniformly. A criteria for 
the choice of h\ was a minimum of L 2 -norm of the 
residual in fl. 

4. Numerical Experiments for 
a ID Eigenvalue Problem 

Accurate approximation of eigenvalue problems is 
essential for bifurcation analysis of PDEs. We 
have not found references in literature on the MQ- 
solution of eigenvalue problems. We therefore 
present here results for an eigenvalue problem for 
ID Laplace operator. For details on the MQ dis- 
cretization see Sec. 3. This is a scalar problem 

-u" = \u, u(0) = u(l) = 0 , (12) 

that has the exact solution: 

(A m , U m (x )) = ((7rm) 2 , sin(Trmx)) , 

m = 1,2,... 

where (A m , U m {x )) is the mth eigenpair of (12). 
Introduce the mesh x n = nh, n = 0, 1 , N, 
h _ j/jv, and consider the standard second-order 
finite difference (FDM) discretization of (12): 

Un-t-l — 2u n + U n _i _ . 

h? ' (13) 

n = 1, . . . , N - 1 , u 0 = u/v=0. 


The corresponding approximate eigenpairs are 
given by 




7T771 

sin lv 

\ 

(>£,. O = 

, , irm 

4iV sin w , 

7r2m 

sin N 




-n(N - l)m 




L sm N 

J 


771 = 1, . . . , N — 1 . 

We also solved (12) using the MQ discretization for 
several values of the number K of internal nodes. 
Denote by (A"«, m = L ■ • • , K the corre- 

sponding approximate eigenpairs. 

The results of our computations are summa- 
rized in Table 1. We use the notation cf Q . s\ for 
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the relative errors in A^, A^, respectively, and the 
notation for the Loo-norm error in U% Q . For 
each MQ solution we provide a comparison with the 
FDM solution that has a sufficient number of nodes 
to give the same accuracy for A i as the MQ method. 
In Part (a) of the table we use the uniform node 
distribution for the MQ method. Part (b) of the 
table shows that the accuracy of the MQ method 
can -be significantly improved by adapting the node 
distribution: we moved only two nodes adjacent to 
boundary to reduce their distance from the bound- 
ary to h x = h/4 (while the remaining nodes are 
distributed uniformly). 

One can see that the MQ method can give a 
highly accurate solution with a small number of un- 
knowns, 10-100 times smaller than the number of 
unknowns in the FDM for the same accuracy. 


5. Numerical Experiments for 
ID and 2D Elliptic PDEs 

We present several examples of continuation of so- 
lutions to systems of nonlinear ID and 2D ellip- 
tic PDEs. Each problem is discretized by the MQ 
method described in Sec. 3. We then perform con- 
tinuation of the resulting system of algebraic equa- 
tions (11) with AUTO 97. The principal goal of our 
examples is to assess the accuracy of the detection 
of bifurcation points. We compare our results with 
some published results and, in one case, the results 
of our computations with an example in AUT097 
and CONTENT. We will use throughout the nota- 
tion K for the number of unknowns in a particular 
method. For our MQ method K = nxN, where n is 
the dimension of the system and N is the number of 
internal nodes. We denote by MQ(u) and MQ(nu) 
our MQ method with the uniform and nonuniform 
node distributions, respectively. 

Example 1. ID Gelfand-Bratu equation. This is 
a scalar problem 

U " + Ae u = 0, in f! = (0, 1) , 

(14) 

u(0) = u(l) = 0 , 

that appears in combustion theory and is used as 
the demo example exp in AUT097 [Doedel et al., 
1997] (fifth-order adaptive orthogonal spline collo- 
cation method) and demo example brg in 
TENT [Kuznetsov & Levitin, 1995-1997] (third- 
order adaptive finite difference method). There is 
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Table 1. A ID eigenvalue problem: comparison of results for eigenvalues, results for eigenfunctions. 


(a) The MQ method with a uniform node distribution for K = 5, 7 and 9. 


m 

A m (Exact) 

'Vn Q . K = 5 

Rel. Err. 

Rei. Err. e^J Q 

Rel. Err. e x , K = 

47 

1 

9.86961 

9.86596 

3.7 x 10~ 4 

3.7 x 10" 4 

3.7 x 10 -4 


2 

39.4784 

39.6492 

4.3 x 10 -3 

5.2 x 10 -3 

1.5 x 10” 3 


m 

A m (Exact) 

Am Q , A = 7 

Rei. Err. 

Rel. Err. £y Q 

Rel. Err. e h x , K = 

: 76 

1 

9.86961 

9.86821 

1.4 x 10 -4 ' • 

■ 9.9 x 10 -5 

1.4 x 10 -4 


2 

39.4784 

39.4738 

1.2 x 10~ 4 

1.8 x 10~* 

5.7 x 10" 4 


3 

88.8264 

89.3648 

6.0 x 10 -3 

1.1 x 10 -2 

1.3 x 10" 3 


m 

A m (Exact) 

A“ Q , K = 9 

Rel. Err. 

Rel. Err. 

Rel. Err. e$J t K = 

117 

1 

9.86961 

9.86901 

6.0 x 10 -5 

5.0 x 10" s 

6.0 x 10“ 5 


2 

39.4784 

39.4846 

1.6 x 10 -4 

2.1 x 10 -4 

2.4 x 10 -4 


3 

88.8264 

89.1667 

3.8 x 10" 3 

7.3 x 10 -3 

5.4 x 10 -4 


4 

157.913 

159.689 

1.1 x KT 2 

2.5 x 10 -2 

9.6 x 10" 4 



(b) The MQ method with nonuniform node distribution for K 

= 7 and 9. 


m 

A m (Exact) 

\Z! q ,K = 7 

Rel. Err. e x Q 

Rel. Err. e™' Q 

Rel. Err. s x , K = 

3477 

1 

9.86961 

9.86961 

6.8 x 10~ 8 

3.0 x 10“ 6 

6.8 x 10“ 8 


2 

39.4784 

39.4782 

3.2 x 10 -6 

3.0 x 10~ 4 

2.7 x 10~ 7 


3 

88.8264 

88.8139 

1.4 x 10 -4 

6.5 x 10 -4 

6.1 x 10~ 7 


m 

A m (Exact) 

. K = 9 

Rel. Err. ef Q 

Rel. Err. t% Q 

Rel. Err. e x , K = 

950 

1 

9.86961 

9.86960 

9.1 x 10 -7 

2.3 x 10~ 6 

9.1 x 10~ 7 


2 

39.4784 

39 4733 

1.4 x 10" 6 

2.0 x 10 -5 

3.6 x 10~ 6 


3 

88.8264 

88.8241 

2.6 x 10 -5 

1.8 x 10 -4 

8.2 x 10" 6 


4 

157.913 

157.882 

1.9 x 10~ 4 

1.8 x 10~ 3 

1.5 x 10" 5 



a limit (fold) point on the solution curve. We take 
the value of A at the limit point found from demo 
exp ( K > 50) as exact. Table 2 shows the compar- 
ison between numerical results in [Davidson, 1997], 
our numerical results and our experiments with 
CONTENT. 


used in the literature as a standard model for bi- 
furcation analysis, see e.g. [Schaeffer &: Golubitsky, 
1981; Golubitsky & Schaeffer, 1985; Dangelmayr, 
1987; Chien et al., 1997, Eq. (24)] and [Mei, 1997], 
A stationary bifurcation occurs [Chien et al., 1997, 
Eq. (24)] at 


Example 2. ID Brusselator problem, a well- 
known model system for autocatalvtic chemical re- 
actions with diffusion: 

u" — {b + l)u 4- u 2 v + a = 0 , 

^■v" -i- bu - u 2 v = 0, in D = (0, 1) , (15) 

l 2 

u(0) = u(l) = a, v(0) = v(l) — - , 

proposed in [Lefever & Prigogine. 1968]. This prob- 
lem exhibits a rich bifurcation scenario and has been 


b n = 1 + j^a 2 + 

d 2 


7 r 2 n 2 


d\ 4- 


7 r 2 n 2 d 2 


> 0 . 


For l = di = 1, d 2 = 2, a = 4, n = 1, 2 this 
gives simple bifurcations: b\ = 9 4- tt 2 ■+■ 8/7r“ = 
19.680174, b 2 = 9 -I- 4tt 2 4- 2/tv 2 = 48.681060, 
correspondingly. For the second-order central dif- 
ference method with uniform mesh of 41 mesh 
points ( K = 80 unknowns), the corresponding ap- 
proximate bifurcation points were found in [Chien 
et al, 1997, Sec. 6.1]. Table 3 shows the comparison 
between analytical, numerical results [Chien et al., 
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Table 2. 

ID Gelfand-Bratu equation: The limit point comparison. 



(a) Results for the MQ method correspond to a uniform node distribution. 



[Doedel et al., 1997], 
Exact 

[Davidson, 1997], 
K = 800 

MQ(u), K = 5 

MQ(u), K = 7 

MQ(u), K = 9 

A 

rel. error 

3.513831 

3.5137 
3.7 x IO -5 

3.512609 
3.5 x 10 -4 

3.514224 
-1.1 x IO' 4 

3.514047 
-6.1 x 10' s 

(M R pstilts for the MQ method correspond to a nonuniform node distribution. 


[Kuznetsov & 
Levitin, 1995-1997], 
K = 50 

[Kuznetsov 
Levitin, 1995-1997], 
K = 500 

MQ(nu), K = 5 

MQ(nu), K = 7 

MQ(nu), K -9 

\ 

rel. error 

3.51145 
6.8 x 1CT 4 

3.51380 
8.8 x IO' 6 

3.514010 
-5.1 x 10 _s 

3.513809 
6. 3 x 10' 6 

3.513828 
8.5 x 10~ 7 


61 

rel. error 


b 2 

rel. error 


Table 3. ID Brusselator equation. Comparison for the bifurcation points. 

(a) Results for the bifurcation point hi- 


Exact 


19.680174 


Exact 


48 681060 


[Chien et al., 
1997], K - 80 


MQ(u), K = 10 MQ(u), K = 14 MQ(u), K = 18 


19.67547 
2.4 x 10' 4 


19.67366 
3.3 x 10~ 4 


19.67786 
1.2 x 10' 4 


19.67919 
5.0 x 10' 5 


(b) Results for the bifurcation point fo- 


[Chien et al., 
1997], K = 80 


MQ(u), K = 10 MQ(u), K - 14 MQ(u), K 18 


48. 6004 
1.7 x IO' 3 


48.57476 
2.2 x lO" 3 


48.63168 
1.0 x IO -3 


48.65605 
5. 1 x 10~ 4 


1997, Sec. 6.1] and our numerical results for values 
of bi and b 2 at simple bifurcation points. 

Example 3. Pattern formation in a ID system 
with mixed boundary conditions [Dillon et al., 
1994]. 

A-u" 4 - 3 — ku — uv 2 - 0 , 
ujI 2 

s — v" 4 - KU - uv 2 - v = 0, in 9 = (0, 1) 

l 2 ( 16 ) 

B\ ^ = p(l - dl){&3U S - u ) - 
on 

S0, — = 6p(l - 9 2 ){9 2 v s - v), on 89 = {0, 1} . 

“ On 

Here Q, e iO. 1], i = 1. 2, 3. are homotopy parame- 
ters. For di = IO' 5 . ^ = IO" 2 . & = °- 14 - 3 = 1.0, 


k = 0.001, (01, 0 2 , 9 3 ) = (1, 1, o) (Neumann prob- 
lem). Equation (16) was discretized by the second- 
order central difference method with equidistant 
mesh of 41 mesh points {K = 80 unknowns). Ta- 
ble 4 [Dillon et al., 1994, Table 1] shows a com- 
parison between analytic and numerical results for 
values of Z at simple bifurcation points. 

Our numerical results (MQ(nu) method) with 
X — 18 coincide with the analytic results above. 

Example 4. 2D Gelfand-Bratu problem 

Aix 4- Ae u = 0, fi = (0, 1) x (0, 1), 
u|an = 0 . 

This problem was studied in [Schwetlick et al.. 1996] 
and [Doedel & Sharifi, to appear]. In [Schwetlick 
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Table 4. A ID pattern formation problem, comparison for simple bifurcation points. 


[[Dillon et al., 1994], numerical] 
[[Dillon et al., 1994], analytic] 
MQ(nu) 


0.047 

O.OSO 

0.093 

0.0465 

0.0793 

0.093 

0.0465 

0.0793 

0.093 


0.159 

0.140 

0.238 

0.159 

0.140 

0.238 

0.159 

0.140 

0.238 


0.186 

0.317 

0.232 

0.186 

0.317 

0.233 

0.186 

0.317 

0.233 


Table 5. 2D Bratu equation, results for the limit point. 


(a) Uniform node distribution. 



[Schwetlick et al., 1996], 
225 < K < 3025 

MQ(u), K = 25 

MQ(u), K = 49 

MQ(u), K = 81 

A 

rel. error 

not reported 

6.8359 
4.1 x 10 -3 

6.8183 
1.5 x 10 -3 

6.8119 
5.6 x 10 -4 

(b) nonuniform node distribution. 


[Doedel & Sharifi, to appear], 
Exact 

MQ(nu), K = 25 

MQ(nu), K = 49 

MQ(nu), K = 81 

A 

rel. error 

6.808124423 

6.793248 
-2.1 x 10 -3 

6.807978 
-2.2 x 10 -5 

6.808232 
-1.6 x 10 -5 


et al., 1996] it was discretized with the second- 
order central difference method with uniform mesh 
and then continued using Implicit Block Elimina- 
tion based on Recursive Projections. A limit point 
was detected for some value of A (not reported in 
the paper), and spurious limit points were detected 
for K = 961, 1521, 2209, 3025 and A sufficiently 
small. In [Doedel k Sharifi, to appear] the prob- 
lem was discretized with a high-order orthogonal 
spline collocation method with sparse Jacobian. We 
reproduced the bifurcation diagram in [Schwetlick 
et al., 1996]. Table 5 gives the values of A at the 
limit point computed by the MQ method. The ex- 
act location of the limit point is assumed to be at 
the value of A obtained in [Doedel & Sharifi, to ap- 
pear] on a 16 x 16 mesh with 4x4 collocation points. 
See Sec. 6 for a discussion of the operation count. 

Example 5. 2D Brusselator problem. 

^ Au - (b + l)u u 2 v + a = 0 , 

l 2 

^ Av + bu - u 2 v = 0, in Q = (0, 1) x (0. 1) , 

, b 

u |an= a, v |aa= - . 

18 ) 


A stationary bifurcation occurs [Chien k Chen. 
1998, Eq. (2.26)] for 

= 1 + 1" 2 + ^ (t + " 2 ) 


7r 2 d2 \m 2 + l 2 n 2 J 

For l = dx = 1, d 2 = 2 , o = 4 , (m, n) = 
(1, 1), (m, n) = (2, 2) this gives simple bifurca- 
tions: 6i,i = 9 + 2 tt 2 +4/tt 2 , 6 2 ,2 = 9 + 8 tt 2 4- l / it 2 , 
correspondingly. For the second-order central dif- 
ference method with equidistant mesh of 21 mesh 
points, the corresponding approximate bifurcation 
points are found in [Chien & Chen, 1998, Sec. 5]. 
Tables 6 and 7 show comparisons between ana- 
lytical, numerical results [Chien k Chen, 1998. 
Eq. (2.26)] and our numerical results for values of 
611 and 62,2 a t simple bifurcation points. 

A Hopf bifurcation occurs [Chien k Chen, 1998, 
Eq. (2.26)] for 

b m>n = 1 + a 2 + (dr + d 2 ) + n 2 ^j tt 2 

for some m, n, and l large enough. For l = 10. 
di = d 2 = 1, a = 10, (m, n) = (1, 2), this gives a 
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Table 6. 2D Brusselator equation: comparison for the bifurcation points, a uniform node distribution for MQ 
method. 

(a) Results for the bifurcation point b\. 


b\, i 

rel. error 


^> 2 . 2 

rel. error 


Exact 


29.144494 


88.058156 


[Chien & Chen, 1998], 

K = 800 MQ(u), K = 50 MQ(u), K = 72 MQ(u), Kj-98 


29.104774 
1.4 x 10 -3 


29.16280 
-6.3 x 10 -4 


29.17050 
-8.9 x 10“ 


(b) Results for the bifurcation point fa. 


87.47325 
6.6 x 10“ 3 


87.61578 
5.0 x 10 -3 


87.86924 
2.1 x 10 -3 


29.16062 
-5.5 x 10 -4 


[Chien & Chen, 1998], 

Exact K = 800 MQ(u), if = 50 MQ(u), if = 72 MQ(u), if - 98 

88.00143 


6.4 x 10' 


Table 7. 2D Brusselator equation: Comparison for the bifurcation points, a nonuniform node 
distribution for the MQ method. 

(a) Results for the bifurcation point 6i, 1 ■ 


61.1 

rel. error 


b 2, 2 

rel. error 


Exact MQ(nu), if = 50 MQ(nu), if = 72 MQ(nu),if=98 


29.144494 


29.14621 
-5.9 x 10 -5 


29.14726 
-9.5 x 10 _s 


(b) Results for the bifurcation point 61,2. 


88.058156 


88.15470 
-1.1 x 10" 


87.93391 
1.4 x 10" 3 


29.14431 
6.3 x 10" 6 


Exact MQ(nu), K = 50 MQ(nu), K = 72 MQ(nu), if = 98 


88.07288 
-1.7 x 10" 4 


Table 8. 2D Brusselator equation, results for the Hopf bifurcation point. 



Exact 

MQ(u), K - 50 

MQ(nu), K = 50 

MQ(u), K = 72 

MQ(u), K =98 

b\,2 

rel. error 

180. 15 

181.8625 
-9.5 x 10" 3 

180.7880 
-3.5 x 10" 3 

181.0696 
-5.1 x 10" 3 

180.492 
-1.9 x 10" 3 


Hopf bifurcation at b h2 = 101 ^ 2 ( ( 1 / 100 ) H- 2 2 )tt 2 = 
180.15. see Table 8. 


6. Discussion 

1. We have presented the results of our exper- 
iments with the continuation of solutions to ID 
and 2D nonlinear elliptic PDEs discretized by the 
MO method. We use a small number of unknowns 
and achieve a high accuracy for detected bifurca- 
tion points in our examples. Here are some sample 
results. 


(i) For the limit point in the ID Gelfand-Bratu 
equation, the MQ method with nine unknowns 
gives the relative errors 6.1 x 10“° and 8.5 x 
l(p 7 for the uniform and nonuniform node dis- 
tributions, respectively. The relative error in 
the third-order finite difference method with 
500 nodes is 8.8 x 10~ 6 . 

(ii) For the two bifurcation points in the 2D Brus- 

selator problem, the MQ method with 98 un- 
knowns gives the relative errors 5.5 x 10 , 

6.4 x 10~ 4 for the uniform node distribution 
and 6.3 x 10” 6 , 1.7 x 10” 5 for the nonuniform 
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node distribution. The corresponding rela- 
tive errors in the second-order finite differ- 
ence method with 800 nodes are 1.4 x 10 3 . 
6.6 x 10“ 3 . 

(iii) For the first eigenvalue in the eigenvalue prob- 
lem for the ID Laplace operator with nine un- 
knowns gives the relative error 6 x 10 and 
9 x 10 -7 for the uniform and nonuniform node 
distributions, respectively. This is equivalent 
in accuracy to 117 and 950 node solutions, re- 
spectively by the second-order finite difference 
method. 

2. As we noted in Introduction, the MQ method 

leads to systems with full matrices. Solving a re- 
lated linear system for the number of nodes M x M 
in 2D with a full M 2 x M 2 matrix by Gaussian 
elimination takes 2/3M^ + 0(i\/ 4 ) operations. By 
comparison, a band solver would take 0(M 4 ) op- 
erations, and a collocation method on a square 
[Doedel. 1998; Doedel &: Sharifi, to appear] would 
take 62 p 3 M 3 , where p is the number of matching 

points at an edge of a finite element [Doedel, 1998]. 
Further work is required to carefully compare the 
costs of solving linear systems and the correspond- 
ing eigenvalue problems arising in discretizing el- 
liptic PDEs by the MQ method and by the finite 
difference, finite element, and collocation methods. 

3. An increase of the number of unknowns and 
especially the shape parameter result in a better 
accuracy but also in a larger condition number of 
the operator T mapping the nodal values of the so- 
lution onto the expansion coefficients. This condi- 
tion number is a limiting factor in our experiments. 
In fact, to reach high accuracy for the limit point 
in the 2D Gelfand-Bratu problem (e.g. the rela- 
tive error 1.6 x 10 -5 with 81 unknowns), we had 
to use quadruple precision. This is a temporary 
fix. as it considerably slows down computations. In 
future, we plan to implement the ideas of Kansa 
et al. [1990b] to circumvent this ill-conditioning 
problem. 

We also found that even a simple adaptation 
of the nodes adjacent to the boundary can lead to 
a dramatic improvement of the accuracy in detect- 
ing bifurcation points. Adaptive choice of the shape 
parameter is another way to improve the accuracy 
that we plan to investigate. 
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Abstract 


Convective transport of heat and constituent components 
is dominated by buoyancy driven convection in many ter- 
restrial crystal growth situations. The character of natural 
buoyant convection in non-uniformly heated, rigidly con- 
tained inhomogeneous fluids can be drastically altered by 
vibration of the container. Therefore, vibrational induced 
flow can potentially be used to influence and even control 
transport in some crystal growth situations. 

A parametric numerical investigation of 3D thermovibra- 
tional buoyancy-driven flow in differentially heated cylin- 
drical containers has been conducted to investigate ther- 
movibrational transport regimes in Bridgman-type sys- 
tems. The objective of the work is to assess the feasibility 
of the use of vibration to suppress, or control, convection 
in order to achieve transport control during cry stal growth. 

The formulation of a model for this problem is outlined, 
numerical method is described and its application to the 
study of investigation of thermal vibrational flows is dis- 
cussed. Two types of vibration are considered: transla- 
tional. and circularly polarized. The results tor flows in- 
duced bv g-jitter and selected results for the cases of lon- 
gitudinal and lateral vibrations are presented. 


"Copvrighi 02000 by A. I. Fedoseyev. Published b> ihe American 
Insutuie of Aeronautics and Astronautics. Inc., with per.-.ission. 


1 Introduction 

It is generally recognized that oscillator)', or pulsative. 
flow significantly alters the transfer of mass, heat and mo- 
mentum in fluid systems. For certain experiments and 
operating conditions, vibrations are expected to ha\e a 
significant influence on heat and mass transfer onboard 
the International Space Station (see for example, the re- 
cent ESTEC Workshop proceedings [1])- Available flight 
experiment data clearly show that, once initiated by g- 
ji U ef\ the effects of convective flows can persist for long 
times even when the g-jitter disturbance (and consequent 
flow) were short-lived [2]-[7]. 

Control of convective transport continues to be an impor- 
tant aspect of crystal growth research. Several groups are 
actively pursuing control of convection using static and 
rotating magnetic fields. However, magnetic fields cannot 
be used for flow control in melts and solutions that are 
poor conductors. Flow control through vibration or vibro- 
convective mixing may offer an attractive alternative in 
such cases. 

Recent works have shown that the character of natural 
buoyant convection in non-uniformly heated, rigidly con- 
tained inhomogeneous fluids can be drastically altered by 
vibration of the container. A review and relevant theoret- 
ical and experimental research can be found in publica- 
tions [1]-[13]. Thus, vibrational induced How can poten- 
tially be used to influence and even control transport in 
some crystal growth situations. A practical quantitati\e 
understanding of vibrational convection as a control pa- 
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rameter in crystal growth situations is currently not avail- 
able. The objective of the work is to assess the feasibility 
of the use of vibration to suppress, or control, convection 
in order to achieve transport control during crystal growth. 

2 Problem formulation and numeri- 
cal method 

Buoyancy driven vibro-convective motion occurs when 
oscillatory displacement of a container wall induces the 
acceleration of a container wall relative to the inner fluid. 

If the fluid density is nonuniform, fluid motion may ensue. 
The magnitude of this motion depends on the orientation 
of the vibrational direction with respect to the local den- 
sity gradients. It should be noted that even in case of a 
constant density fluid subject to spatially nonuniform vi- 
bration. fluid motion can also occur (tor example, angular 
vibration [ 11 ]). 

To properly investigate influence of translational and cir- 
cularly polarized vibration necessitates the use of the full 
3D equations governing the transport of heat, mass and 
momentum. Selected examples ot our ongoing work on 
this topic are outlined below. 

We consider a purely thermo-vibrational convection in a 
differentially heated cylindrical cavity. The fluid is taken 
to be Newtonian, and the Boussinesq approximation is as- 
sumed to hold. The calculations were performed tor iden- 
tification and characterization of thermovibrational flow 
and are part of an ongoing project involving flow visual- 
ization model experiments being conducted by Feigelson 
[ 101 . 

2.1 Governing equations 

Translational vibration corresponds to a linear displace- 
ment such as, tor example, u — dcoscor, where d is a 
real vector giving the displacement magnitude and to is 
the frequency. In this case the ampoule is displaced back 
and forth upon the same line. Polarized vibrations are 
characterized by a displacement u = /^{dc' 11 *}, where 
d = di - fdi. Here the instantaneous vibration direction 
rotates in the polarization plane defined by the real vec- 
tors di and d 2 . A sketch showing both translational and 
circular polarized vibrations is presented in Fig. 1 . In a ref- 
erence frame fixed to a vibrating ampoule, these types of 
vibrations result in the following form of the momentum 
equation: 


+ (wy V = -Wp + PrV 2 V + Ra T Pr • (O + aC)fi*+ 

Ra* r Pr • (0 + aC)f(£2, r) (1) 

while the continuity equation: 

V • V = 0 (2) 

energy equation: 

^ + (VV)G = V 2 © (3) 

dt 

and species transport equation: 

— + (VV)C = PrSc~ l V 2 C (4) 

dt 

where length, time and velocity are scaled, respectively, 
by Ro , Rq/k and k//?o- Here Ro is the ampoule radius and 
k is the thermal diffusivity. The nondimensional concen- 
tration and temperature, are given by 0 and C, respec- 
tively. The function f(n,f) is the acceleration of the vi- 
brating ampoule and £2 = co/?q/k is a dimensionless fre- 
quency; n s = (sin<t>,0,cos<t>). The Prandtl, Schmidt, ther- 
mal and solutal Rayleigh and vibrational Rayleigh num- 
bers and the buoyancy ratio are, respectively, given by 

Pr*=Z,Sc=i,Ra T = ^ 7 ^ , Ra s = aRa T . 

a = p ( .c./(PAr),«4 = ‘^^ 

Here p and p t are the thermal and solutal expansion co- 
efficients and AT, c~, g, to, d, k, v, D are the character- 
istic longitudinal temperature difference, reference con- 
centration in the melt, gravitational acceleration, vibra- 
tional displacement amplitude and frequency, direction of 
gravity, kinematic viscosity and solute diffusivity, respec- 
tively. The dimensionless number Ra* r is the vibrational 
Rayleigh number and Ra * = o.Ra* r . 

The boundary conditions for eq. (l)-(3) are: (i) the non- 
slip condition for the velocity, V = 0, on the walls, and 
(ii) the given wall temperature distribution, 0 = 0 W . An 
undefined constant in the pressure field is excluded by set- 
ting p{xo,yo,zo) = 0 at some location (*o,yo,3>)- 

Equations (l) are solved together with the equations gov- 
erning heat and species transfer ( 3 ), (4) and the condition 
that the velocity is divergence free ( 2 ). 
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2.2 Finite element solution method 

The equations arc solved in primitive variable form 
(velocity-pressure, temperature, concentration, etc.) us- 
ing a Finite Element Method. 

Finite element method with regularization for the 
Navier-Stokes equations 

This method (RNS) was proposed in [151, [16] high 
Rc number flows. It was shown that such a regularization 
works also well in case of flows with thin boundary' lay- 
ers, even with few mesh nodes placed inside the boundary 
layer [17], [181. For the considered problems, the conti- 
nuity equation (2) is modified as follows 

V- V = xV*(V/>-F) (5) 

where i is a small regularization parameter, and F is a 
body force in eq. (1) for the thermo-vibrational convec- 
tion. For x -> 0 we approach the original equation (2). A 
boundary condition for the pressure on the wall is 

(V/7 — F) ■ n = 0, (6) 

where n is a unit wall normal vector. Eq. (5) and (6) 
present the main feature of this method, and ensure the 
balance of the component of the force normal to the region 
boundary. 

Another advantage is that this approach allows to use the 
same order finite element approximation for the veloc- 
itv. pressure, temperature and concentration with all un- 
knowns located at the same nodal points. For a justifica- 
tion of this regularization one can be referenced to the new 
hydrodvnamic equations proposed in [191 that have simi- 
lar fluctuation terms. Ldhner has also shown that similar 
terms actually appear in the discrete equations as a result 
of different order finite element approximations used for 
interpolation of velocity and pressure [20]. 

The continuity equation (5) and momentum equations (I) 
are solved simultaneously at each time step. This elimi- 
nates manv problems related to boundary conditions and 
places only slight limitations on the time step size for tran- 
sient problems (due to the physical nature of the problem). 

3D CFD software 

We implemented the above 3D model of vibro-convective 
buovancy-driven flow in differentially in the FEMINA/3D 


code [14], The regularization proposed makes a solution 
procedure very efficient. 

A high accuracy solution iterative CGS-like method using 
preconditioning by high order incomplete decomposition 
has been implemented. This allowed us to obtain high- 
precision solutions with accuracy up to 1CT 9 . The precon- 
ditioning also reduced the computation time by one to two 
orders of magnitude and the memory size by a factor of 8 
for3D flows compared to currently available commercial 
codes (e.g. CFD2000) [24]. This approach makes it possi- 
ble to solve large time-dependent problems (up to 300.000 
unknowns) with reasonable computation times. A typical 
solution time for a transient problem is few hours on a SGI 
02 workstation. 

2.3 Benchmarks 

This code was carefully tested on benchmarked exper- 
imental, theoretical and numerical data for a variety 
of 2D/3D viscous and thermo-convective flow problems 
[18], [16], [17]. Here we present some selected examples: 

Three dimensional thermal convection in a cylinder 

The method bees applied to the problem of convective 3D 
flow in differentially heated horizontal cylinder. The ex- 
perimental data by Bogatirev et al [21] have been used 
for comparison. These data have been obtained during 
ground tests for the device, the thermal convection sen- 
sor, before it was flown on Mir station. Numerical re- 
sults from the 3D finite volume simulations by Bessono\ 
[22] have been also used for comparison. The temperature 
distribution on a cylinder wall was (i) linear temperature 
profile, and (ii) computed using a real, finite wall con- 
ductivity (adjoint problem). The body force in eq. (1) 
is f = (0,0 ,Ra T PrO), no vibration was applied. The 
Rayleigh numbers is in the range from 10*’ to 1.2* 10\ 
The value of x* used was 10“ 7 to 10~\ and it did not 
change noticeably the results. Results are shown in Fig. 

An agreement with the experimental data for Raj >4-10 
is quite good in the case (ii), when a real finite wall con- 
ductivity has been taken into account. 

Two dimensional and three dimensional lid-driven 
cavity problem 

We compare our results with experimental data obtained 
by Koseff and Street [23] for isothermal flow at Re = 
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3200 and 10 4 . We solved equations (1) and (5) numeri- 
cally for the unknowns (V, p) in the 2D region {.x,z) = 

[0, 1][0, 1]. The 3D version of the problem corresponds 
to’ the Koseff & Street experiment [23] with the domain 
(.v,.v,c) = [0, 1][0, 3][0, 1]. The boundary conditions are: 

V = (»,v,w) = (1,0,0) at the driven lid (z= 1); V = 0 on 
the other walls, and, (Vp • n) = 0 on all the walls. The 
undefined constant in the pressure field is eliminated by 
setting pressure p = 0 at (x,y,c) = (0,0,0). 

Fia. 3 and 4 present the experimental measurements and 
numerical results for u and w velocity components ob- 
tained by our method. RNS, by solving the original NS 
equations (1), (2), and by using a 2D k — E model with 
commercial code [24]. The experimental data are shown 
for the symmetry plane (x,z) at y=1.5 along the lines x 
=0.5 and z =0.5 of the 3D cavity, that has relative dimen- 
sions (x:y:z)= 1:3:1. The experimental points correspond 
to time averaged values of the measured velocities [23]. 

The results obtained with our model, eq. (1). (5) and 
(6) are in good agreement with the experimental data for 
Re=3200 and 10,000 and are an improvement upon previ- 
ous results obtained by solving the NS equations (1), (2). 

Magnetic field suppression of convective flow 

Numerical solution for thermal convection flows in a 
semiconductor melt with strong static magnetic field ap- 
plied is presented in Fig 5. Although the generated flows 
have extremely low velocity because ot the large Hart- 
mann numbers (Ha = 20 to 2000). the numerical solu- 
tion of the governing equations involved is very compli- 
cated due to the thin boundary layers. Different numerical 
methods have been tested for the solution of this problem 
[17]. The best results have been obtained with the pre- 
sented RNS approach. It can provide the numerical solu- 
tions in a w'ide range of Ha numbers (up to 10 4 ), while 
other methods failed for Ha > 20. The RNS results com- 
pare favorably with the asymptotic theoretical solutions. 
Fig. 5c. 

3 Thermal vibrational convection. 
Results and discussion 

A parametric study of translational and circular polarized 
vibrations under typical microgravity and terrestrial con- 
ditions for typical semiconductor melts was performed. A 
snapshot of a typical flow pattern tor translational vibra- 
tion is presented in Fig. 6. Even in the total absence ot 


gravity the vibrations have resulted in detectable flows. 

For the cases examined, the temperature distribution re- 
mains almost unperturbed (due to the low Pr and weak 
flow strength). 

The anele between the direction of vibration and the am- 
poule has been studied for translational vibrations in the 
presence of an axial temperature gradient. At high fre- 
quencies and when the angle is zero, no influence of the 
vibratipn on the flow was observed, even when vibrational 
the Raleigh number is very high. The maximum observed 
effect corresponds to an angle of 90 degrees. Here trans- 
port is significantly enhanced. 

The influence of vibrations on heat and mass transfer be- 
comes significant for oxide melts due to their low thermal 
diffusivity (Pr ~ 10). These flow patterns are shown in 
Fie. 7 for the case of circular polarized vibration. Initially 
(at" time r = 0), the species concentration was c = 1 at the 
lower quarter of the cylinder and c = 0 elsewhere. The 
evolution of the species concentration (process of mixing) 
is shown in Fig. 8 together with minimum and maximum 
values of velocity (for the whole domain) components. 
Complete mixing occurs in about ten seconds. The heat 
transfer (local Nusselt number at the top and the bottom) 
is also enhanced by about an order of magnitude. If the 
frequency of vibration is high, of the order of 100/7: (for 
fixed Ra ), then the changes in heat and mass transfer due 
to vibrations become less significant. This corresponds to 
earlier experimental observations [7], [8]. 

Our results show that both translational and circular po- 
larized vibrations can cause average melt flow for a range 
of parameters typical of practical semiconductor growth. 
For a given vibration amplitude and frequency, circular 
polarized vibrations result in more intensive melt flows 
than translational ones. 

The influence of forced vibration on g-jitter induced flows 
using typical SAMS micro acceleration data from the 
USML-2 mission was also investigated. Motivated by the 
predictions of the averaged equation theory presented in 
Ref. [11], translational vibration was applied parallel to 
the ampoule axis (and thus, the temperature gradient) in 
an attempt to damp unwanted irregular time-dependent 
flow caused by g-jitter Fig. 9a). While the flow varia- 
tion with time becomes more regular, we did not succeed 
in completely suppressing the g-jitter flow (Fig.9b,c). 

We found that the use of the same amplitude vibration 
in the direction orthogonal to the ampoule axis is more 
effective. This induces intensive thermal vibrational flows 
and flow disturbances due to g-jitter become practically 
insignificant (Fig. 9d). 


4 

American Institute of Aeronautics and Astronautics 


Conclusions [5] 

The influence of translational and circularly polarized vi- 
bration in analysis in a model Bridgman melt growth con- 
figuration was investigated. The nature ot the flows pro- ^ 
duced by the types of vibration under consideration ne- 
cessitated the use of the full 3D equations governing the 
transport of heat, mass and momentum. The governing 
equations were solved numerically. Flow patterns for [7] 
translational and circular polarized vibrations and g-jitter 
microaccelerations were analyzed. For translational vi- 
bration. thermovibrational flow is strongly dependent on 
the angle between the vibration direction and the tempera- 
ture gradient. Circular polarized and rotational vibrations 
result in more intensive melt flows than translational ones. [8] 
The simultaneous action of g-jitter and translational vibra- 
tions is currently being study from the viewpoint of using 
applied vibration as a means of flow control. 
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Figure 1: Translational vibration (a), di or d: = 0, and 
polarized vibration (b), di , di = 0; 0 is the angle between 
gravity vector and axis of ampoule, a is the angle between 
direction of vibrations and axis of ampoule 
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Figure 2: Thermal convection in 3D differentially heated 
horizontal cylinder: (a) temperature field in the central X- 
Z plane. Ra = 10 s ; (b) Comparison of temperature differ- 
ence AT versus Ra number with experiment [21]: squares 
- experimental data, solid line - numerical results for per- 
fect wall conductivity, dashed line - results for finite wall 
conductivity; AT is a measured temperature between lo- 
cations 1 and 2 in (a). 


Figure 3: Driven cavity problem. Re = 3200. Comparison 
of horizontal velocity profiles (1-4) for numerical (solid 
and dashed lines) and experimental (squares) results and 
vertical velocity profiles (5-8) for numerical (solid and 
dashed lines) and experimental (triangles) results : 1-NS. 
2-k-e model. 3-RNS (2D), 4-RNS (3D); 5-RNS (3D), 
6-RNS (2D), 7 - k - e model, 8-NS. 



Figure 4: Driven cavity problem, Re = 10,000. Com- 
parison of horizontal velocity profiles (1-3) for numerical 
(dashed lines) and experimental (squares) results and ver- 
tical velocity profiles (4-6) for numerical (dashed lines) 
and experimental (triangles) results: 1-NS, 2 - k - £ 

model, 3-RNS (2D). 4-RNS (2D), 5 - k - £ model, 6-NS 
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THERM*. CONVECTION 



Figure 6: Typical instantaneous 3D melt flow patterns for 
translational vibration at Og, Ra = 0, Ra v = 7.25 * 10 , 
p r = 0.01, co= 1007/c, lateral vibration: velocity compo- 
nents Vx,Vy\Vz. pressure P, temperature T , temperature 
disturbance AT, and velocity module. White color desig- 
nate the maximal value plotted, black one - the minimal 
value. Vibrations are applied along x-directions (horizon- 
taH 



Figure 5: Thermal convection suppression by magnetic 
field. D = 2, H — 1, Ra = 1-25 - 10 5 , Ha = 21 “0 (B = 
S.QTesla): stream function (a), vertical velocity profile 
Vy(x) at y = 0.25 (b), Summary of magnetic field sup- 
presion of the flow for H = 1. D = 2: maximum value of 
horizontal (radial) velocity versus Z? (c). Predicted theo- 
retical asymptotic dependence V mtLX ~ Ha~- is observed 


Figure 7: Instantaneous 3D flow patterns for circular po- 
larized vibration, Ra = 7.25 ■ 10\ Ra v = 7.4- 10 6 , Pr = 15. 
co = 10 Hz : velocity components Vx,Vz> temperature T. 
concentration C, velocity module, and temperature distur- 
bance DT. White color designate the maximal value plot 
ted. black one - the minimal value. 
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Velocity Vz Velocity Vx Concentration 



Figure S: Temporal evolution of (a) species concentration 
C, min. and max. values ot C. (b) velocity extremums Vx, 
and (c) Vz for circular polarized vibration applied to oxide 
melt Ra = 7.25 • 10\ Ra v = 7.4 • 10 6 , Pr = I5,w = 10 Hz 


G-jitler 




Figure 9: Melt flow response to g-jitter (a), g-jitter and 
longitudinal vibration,, Ra v = 7.4* 10 4 , (b) and g-jitte rand 
lateral vibration, Ra v = 7.25 • 10. Min-max values of Vx 
velocity component versus time are presented. 
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Abstract 

A numerical solution for thermal convection flows in a 
semiconductor melt with strong static magnetic field ap- 
plied is presented. Although, the generated flows have ex- 
tremely low velocity, the numerical solution of the govern- 
ing equations involved is very complicated due to the thin 
boundary lavers. Rectangulat cavity with different aspect 
ratios and gravity direction aligned and misaligned with 
the magnetic field vector are considered. Three numerical 
methods are compared. It is shown that the finite element 
approach with regularization can provide the numerical so- 
lutions in a wide range of Ha numbers (up to 10 4 ). The 
results compare favorably with the asymptotic theoretical 
solutions. 


1 Introduction 

The application of magnetic fields is one of the most 
promising approaches for the reduction of convection 
durinc directional solidification of electrical conductive 

-Copyright ©2000 by~A. I. Fedoseyev. Published by the American 
Institute of Aeronautics and Astronautics, Inc., with permission. 


melts (semiconductor crystals). Current technology allows 
the experiments with very strong static fields (up to 80 
KGauss) for which, based on the simple scaling analysis 
in stabilized systems (vertical Bridgman method with ax- 
ial magnetic field), nearly convection free segregation is 
expected, [1]. 

However, the reported experimental studies have yielded 
controversial results [2,3]. The computational methods are, 
therefore, a fundamental tool in the understanding of the 
phenomena accounting during the solidification of semi- 
conductor materials. Moreover, effects like the bending 
of the isomagnetic lines, different aspect ratios and mis- 
alignments between the direction of the gravity and mag- 
netic field vectors can not be easily analyzed with analyti- 
cal methods. 

The reported numerical results are not able to explain the 
experimental data[4,5]. Although the generated flows are 
extremely low, the computational task is complicated be- 
cause of the thin boundary layers [6]. 

Here, three different numerical approaches we have used 
for comparison, : 

(1) The spectral method implemented in [7], 

(2) The finite element method with regularization for 
boundary layers [8], 
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(3) The multiquadric method, a novel method with global 
radial basis functions [9]. 

The results obtained by these three methods are presented 
for a wide range of Hartmann numbers corresponding to 
magnetic fields B from 0.05 to 5.0 Tesla (0.5 to 50.0 
KGauss). Comparison and discussion of accuracy, effi- 
ciency, reliability and agreement with the asymptotic so- 
lution are presented. 


2 Governing equations 


0 = (T -Tcold)/&T,= GL , G is a vertical temperature 
gradient The Prandtl, Rayleigh and Hartmann num- 
bers are, respectively, given by 


Pr= l' Ra =&£-,Ha = UhJl 

K VK V V P 

Here P,g,v,p,aare thermal expansion coefficient, gravita- 
tional acceleration, kinematic viscosity, density and elec- 
trical conductivity, and Bo is the magnetic field intensity, 
e g , eg area unit vectors in the direction of gravitational ac- 
nnd magnetic field. 


The two dimensional steady state thermal convection of in- 
compressible viscous fluid (melt) in a rectangular ampoule 
of diameter D and hight H was considered. More general 
cases are discussed in [17], [18], [19] 

The governing equations (Boussinesq approximation) are : 

The momentum equation: 


The boundary conditions for eq. ( 1 )-(3) are: (i) the non- 
slip condition for the velocity, V = 0, on the walls, and 
(ii) the given wall temperature distribution, 0 = y on side 
walls, and 0b = -0o(jc - L) 2 * at the bottom. The latter 
condition represents a parabolic temperature distribution at 
the bottom boundary. To exclude an undefined constant in 
the pressure field we set p(jc 0 ,J’o) = 0 at arbitrary location 

(*o, >’<))• 


(V V) V — PrV 2 V + Vp = RaPr • 0 ■ e^, + F (1) 
while the continuity equation: 


3 Problem parameters 


V-V = 0 (2) 

and energy equation: 

(VV)0 = v : 0 . (3) 

where length, time and velocity are scaled, respectively, by 
L. L 2 /k and k/L . Here L is the smallest of the ampoule 
diameter D and hight H . and k = ^ is the thermal diffu- 
sivity of the melt. F is a body force due to magnetic field 
(Lorentz force). For an axially symmetric configuration 
the Lorentz force is given by 

F = PrHa : [{\ x e B ) x e B ) = ( PrHa 2 V\ , 0) (4) 

in the two-dimensional case, where V| is the horizontal 
component of the velocity. It does not depend on elec- 
trical potential, because the electrical potential is uniform 
(V<*> = 0), when a vertical magnetic field is applied ( see, 
for example. [1]). This is not valid if the symmetry is 
broken, when the magnetic field direction and the grav- 
ity vector are slightly misaligned. However, to simplify 
the study, we neglect this effect. According to [1], Joule 
heating due to electromagnetic field can be neglected as 
well. The nondimensional temperature 0 is scaled by 


The problem was solved using properties for Germanium 
(Ge) melt exposed to magnetic field having intensity Bo in 
a range from 0 to 5 Tesla (and in few cases up to 50 Tesla); 
the direction of magnetic field vector is axial, e B = (0, 1). 
The corresponding Hartmann number varies from Ha- 0 
to 2170 (and in few cases to Ha = 2.17 • 10 4 ). Ampoule 
geometries considered are (a) L = 1cm , H = 2 cm, and (b) 
L = 2cm, H = 1cm . 

Solutions were obtained for the temperature gradient on a 
side wall G = lOK/cm. The average temperature gradient 
on the bottom (solid/liquid interface) was 0.5 K/ cm. The 
corresponding Rayleigh number is Ra = 1 .24 ■ 10 - . Earth 
gravity g was considered to be (i) aligned to the magnetic 
field direction, g = go(0, - 1) and (i) misaligned by 0.5 de- 
grees. 


4 Numerical methods 

Three numerical methods were used to solve the system of 
eq. ( 1 ) - ( 3 ) for the boundary conditions and parameters 
given above. 
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4.1 Spectral element method 

By using the spectral element method (SEM) code NEK- 
TON [7], the domain was divided into quadrilateral ele- 
ments, refined near the walls (Fig. la), and 8x8 Cheby- 
shev polynomials are used inside each element for field 
variables approximation. Total number of elements in our 
tests was 362, and total number of unknowns was about 
5 • 10 4 . 

4.2 Finite element method with regulariza- 
tion for the Navier-Stokes equations 

The finite element method with regularization for the 
Navier-Stokes equations (FEMR) was proposed in [8] for 
high Re number flows. It was shown that such a regular- 
ization works well in case of flows with thin boundary lay- 
ers, even with few mesh nodes placed inside the boundary 
layer. For the present problem, the continuity equation (2) 
is modified as follows 

V-V = xV-(Vp-F-/?aPr-0-e s ) (5) 

where x is a small regularization parameter. For t-tOwe 
approach the original equation (2). A boundary condition 
for the pressure on the wall is 

( V /7 — F — RaPr • 0 • e g ) • n = 0, (6) 

where n is a unit wall normal vector. Eq. (5) and (6) 
present the main feature of this method, and ensure the 
balance of the component of the force normal to the region 
boundary. 

Another advantage is that this approach allows to use the 
same order finite element approximation for the velocity 
and pressure with all unknowns located at the same nodal 
points. For a justification of this regularization one can be 
referenced to the new hydrodynamic equations proposed 
in [10] that have similar fluctuation term. Similar terms 
have been obtained as a result of the consistent treatment of 
time-advancement for the divergence-equation by Lohner 
(see [11]). Lohner has also shown that similar terms actu- 
ally appear in the discrete equations as a result of different 
order finite element approximations used for interpolation 
of velocity and pressure. 

The numerical solution is insensitive to the value of t. 
The value of x should be chosen within the range 10 8 to 
10 -4 . For a smaller value of x the discrete equations be- 
come nearly incompatible, and numerical solution exhibits 
strong spatial oscillations. 


The simplest linear finite elements were used for numerical 
approximation of velocity, pressure and temperature. We 
used triangular meshes with 40 x 100 and 80 x 100 nodes 
refined near the walls (Fig. lb). The results obtained on 
these meshes are very close, so we used 40 x 100 mesh for 
most of the runs. Total number of nodes and unknowns 
is respectively 4000 and 16 • 10 3 . FEMINA/3D CFD code 
(Finite Element Method IN Applications) [12] was mod- 
ified to implement proposed regularization method. Dis- 
crete finite element equations corresponding to (1), (5), (3) 
were solved together simultaneously by CNSPACK solver 
[12] using the CGS-type iterative technique and high order 
preconditioning by incomplete decomposition. 

4.3 Multiquadric radial basis function 
method 

The Multiquadric Radial Basis Function (MQ) Method is a 
novel meshless collocation method with global basis func- 
tions. The concept of solving partial differential equations 
(PDE) using radial basis functions (RBFs) was introduced 
by Kansa in 1990 [9]. He implemented this approach for 
the solution of hyperbolic, parabolic, and elliptic PDEs us- 
ing the MQ RBFs proposed by Hardy [13], [14] for inter- 
polation of scattered data. 

An RBF is a function that depends only upon the dis- 
tance between a point (.v, y) and a reference node {xj , \ ; ) • 
Among studied RBFs still only the MQ RBFs are proven 
to have an exponential convergence for the function in- 
terpolation [16]. A MQ RBF is given by gj{x, y) = 

/ (x - Xj) 2 + {y - yj) 2 + c 2 } , where cj is called the shape 
parameter. The numerical experiments for parabolic and 
elliptic PDEs by Kansa [9] show high accuracy and effi- 
ciency of the MQ scheme. A brief review on MQ RBF for 
the solution of PDE can be found in [15] and on the RBF- 
PDE Web site [22]. This approach results in modest size 
systems of nonlinear algebraic equations which can be ef- 
ficiently solved by using widely available library routines 
and linear solvers for dense matrices. 

For a given set of N nodes in the domain and at the bound- 
ary, the solution for unknown V, p or 0 is approximated 
as a sum of MQ functions with the coefficients as un- 
known. These coefficients are found by collocating gov- 
erning equations at the internal nodes and boundary condi- 
tions at the boundary nodes. Nonlinear algebraic system is 
solved by Newton method. 

We used up to 25 x 25 uniformly distributed nodes and 
constant shape parameter Cj = co = const for all functions 
. Total number of unknowns is 2500. 
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5 Results and discussion 

5.1 Convection in rectangular cavity with 
H/D=2 

The tests below represent the case of melt zone with aspect 
ratio H /D = 2. 

5.1.1 Flow without magnetic field 

The nondimensional parameters are: Pr = 0.006, Ra = 
1.25 • 10 5 , Ha - 0 D = l,H = 2 and L = D. The temper- 
ature distribution at the bottom is given by 0e = -3.575 • 
10- 3 (l - 4 jc 2 ). The results for the case a = 0 (a is the 
angle between the gravity vector and the vertical axis) are 
shown in Fig. 2. The solution obtained by all three meth- 
ods are close to each other. 

The flow pattern consists of two counter-rotating symmet- 
ric cells, located at lower corners. The perturbation of the 
temperature distribution resulting from 0 = V, set on the 
wall, can not be observed on the plot. Note that the tem- 
perature field is suppressing the flow, which is caused by 
the horizontal temperature gradient at the bottom. 

If the gravity direction is misaligned with the ampoule axis 
by 0.5 deg, flow pattern becomes quite different. A nor- 
mal to the gravity component of the temperature gradient 
becomes a main reason for the thermal convection. A sin- 
gle roll is formed, while the magnitude of melt velocity is 
higher by a factor of two to three. 

5.1.2 Flows under magnetic field 

B = 0.05 Tesla. Ha=21.7. The MQ method did not yeld 
a solution, because the Newton method did not converge 
(since the Jacobian becomes ill-conditioned). 

The solution by the SEM and FEMR methods show no- 
table difference. The SEM solution for the velocity field 


the results still without difficulty, the velocity profiles re- 
main smooth. 

B = 5.0 Tesla, Ha=2170. The results from the SEM com- 
putation showed strong numerical instability. The FEMR 
solution is still quite reasonable: the flow pattern is about 
the same, but flow velocity is about two order of magnitude 
lower than in previous case B = 0.5 Tesla. The boundary 
layers become extremely thin (0.0 1 cm), and therefore al- 
most invisible on a plot (Fig.4). 

In case of a misalignment of gravitational acceleration with 
ampoule axis, the flow pattern changes to one big cell for 
this and all other values of magnetic field considered. 

B = 50 Tesla, Ha=21700. This was done just to test the 
ability of proposed FEMR method, the solution still re- 
mains smooth with even three times more thin boundary 
layer compared to B = 5.0 Tesla. 

Stretching of the stream lines by the magnetic field demon- 
strated in Fig. 5 for all the cases above (aligned gravity 
vector) plus additional case B = 0.005 Tesla, Ha=2.17. 
This effect is mentioned in many papers schematically, but 
computational results were never shown. 


5.1.3 Discussion 

Figure 6 shows the maximum radial velocity calculated, 
using the FEMR method, for different values of the im- 
posed magnetic field B. The maximum of horizontal (ra- 
dial) velocity versus B is presented by few curves, marked 
as “Vr(b)” for FEMR on 40x100 uniform in vertical di- 
rection mesh, by “Vr(/)” for 40x 1 00 refined near all walls 
mesh and by “Vr(d)” for 80x100 mesh refined at the walls. 
Results for misaligned case are presented by the curve, 
labeled as “Vr(a = 0.5)”. We can observe a predicted 
asymptotic dependence ~ Ha~ 2 for all the cases, 
starting from about B = 0.05 Tesla (Ha ~ 20), in accor- 
dance to asymptotic given in [1]. 


exhibit numerical oscillations between the mesh nodes. 
The flow pattern from FEMR is the same as in the ab- 
sence of magnetic field (Fig. 3). Vertical velocity profile 
at y = 0.25 shows a boundary' layer. The flow velocity is 
decreased by about a factor of two. 

B = 0.5 Tesla, Ha=217. The boundary layer becomes very 
thin, and the flow velocity is about two order of magnitude 
lower compared to B — 0. The velocity profiles from the 
SEM computation exhibit spatial oscillation with velocity 
sisn change between mesh nodes. The FEMR can provide 


The main difficulty of this problem is a viscous flow with 
thin boundary layer. Despite the fact that actual flow ve- 
locities are very low and the Reynolds number obtained 
using the computed velocities, is Re ~ 10 _1 to 10" , a big 
value of the Hartmann number results in a relatively small 
coefficient at the highest derivative of the velocity in the 
momentum equation. Solution of such a problem exhibits 
thin boundary layer with the thickness 8 ~ Ha~ [ , and the 
the “equivalent” Reynolds number Re eqv ~ Ha~ , for B=0.5 
Tesla Re eqv = 4.7 • 10 4 and B=50.0 Tesla Re eqv = 4.7 ■ 10 8 . 


4 

American Institute of Aeronautics and Astronautics 



5.2 Thermal convection in cavity with aspect 
ratio H/D=0.5 

The following tests present the case with aspect ratio 
H/D = 0.5, D = 2, H = 1 . The temperature distribution at 
the bottom is given by 0a = -7.150- 10 _3 (1 — a“). Applied 
axial temperature gradient is also G = 70 K/cm. 


the generated flows are extremely low, the computational 
task is very complicated because of the thin boundary layer 
at high Hartmann numbers, Ha » 1. We considered melt 
region geometry with different aspect ratios, and gravity 
direction aligned and misaligned with the magnetic field 
vector. The comparison shows that the finite element ap- 
proach with regularization can obtain stable and reliable 
solutions in a wide range of Ha number, up to 10 . These 
results compare favorably with asymptotic solutions. 


Flow without magnetic field 

The solution obtained by all three methods are also close 
to each other. The flow pattern consists of two counter- 
rotating symmetric ceils, that occupy most of the volume, 
Fig. 7. 

In the case of the gravity misalignment with the ampoule 
axis direction by 0.5 deg, the axial temperature gradient 
becomes a main driving force for the thermal convection. 
This results in the change of flow pattern that becomes con- 
sisting of one big convective cell. 


Flow with magnetic field 


The results are shown in Fig. 8. Again when Ha number 
is high, all the methods except FEMR, exhibit the same 
difficulties as in a case of aspect ratio H/D- 2. A sum- 
mary of the results is shown in Fig. 9. The suppression of 
the flow is essentially same efficient as before with similar 
asymptotic dependences V nuLX ** Ha~~. The velocity pro- 
file in the boundary layer obtained by FEMR is shown in 
Fig. 10. One of the advantages of FEMR is that its so- 
lution remains smooth even at the big change of the slope. 
We can see that the thickness of the vertical boundary layer 
is in agreement with asymptotic solution, 5 ~ Ha~ . The 
tangent velocity derivative at the boundary decreases with 


Ha number as ^ 


' Ha' 


Comparing between Fig. 6 and 9 it is found that misalign- 
ment’s impact on the reducing of the convection is more 
important for aspect ratio 1 . 


Conclusions 


The main difficulty of this problem is that a flow has a very 
thin boundary layer. Despite the fact that actual Reynolds 
number is very low, Re ~ 10“' 1 to 10" 5 6 , a high value of 
the Hartmann number results in a relatively small coef- 
ficient at the velocity Laplacian in the momentum equa- 
tion. Solution of such problem exhibit thin boundary lay- 
ers with related, like for high Reynolds number flows, dif- 
ficulties. That is one of the reasons for the discrepancy in 
the results that numerical studies reported. Both the spec- 
tral method and the multiquadric method with global basis 
functions needs improvement to deal with thin boundary 
layers. Multilevel approximation by Fasshauer [20],[x.l] 
can be one of the ways. 

Numerical solution of these problems by available com- 
mercial CFD codes may be not efficient or not possible. 
Adaptive algorithms can be a promising solution. Devel- 
opment of more accurate and efficient solution methods for 
this problem is necessary. 
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Fieure 1- Mesh used in the spectral element method (a). 362 elements. 8 x 8 Chebyshev polynomials approximation 
inside each element, and mesh used in the finite element method (b). 4000 nodes, 8,000 triangle elements. 
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Figure 2: Thermal convection without magnetic field for geometry l, D - ^ “ 2 * 'hffMR meTho^sTcT 

distribution (a), stream function (b), and vertical velocity profile Vy(x) at y = 0.25 by SEM and FEMR methods (c). 

velocity scale is 0.225 cm/s. 



Figure 3: FED = 2, Ra = 1.25.vl0 5 , Ha = 21.7 (B=0.05 Tesla). Stream functions for gravity vector (a) aligned and (b) 
0.5 degrees misaligned relative to the vertical direction, (c) Nondimensional vertical velocity profile Vy(x) ( . c s 
calculated using the FEMR method 0.25 H from the bottom of the cavity (y/H =0.25) 
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Figure 4’ Thermal convection with magnetic field for geometry 1, D — 1, H — 2, Ra 1.25 10 ,Ha 2170. 0(B 
5.0T«/1): stream function (a), and vertical velocity profile Vy(x) at y = 0.25, velocity profile Vy(x) for misaligned by 

0.5 degree gravity direction at y = 0.5 (c). 
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Figure 5: Searching of .he flow streamlines with increasing magnetic field, (a) B - 0. (b )B = 0.005 Tesla, and (c) .0 
(f) correspond, respectively, to B — 0.05,0.5,5.0 and 50.0 Tesla 

10 

American Institute of Aeronautics and Astronautics 





Figure 6: A summary of magnetic field suppresion of the flow for H = 2, D = 1 - Maximum value of horizontal (radial) 
velocity versus B: Vr(b ) on 40x100 uniform in vertical direction mesh, Vr(f) for 40x100 refined near w s m , 
V r (d ), 80x100 mesh refined at the walls. Predicted asymptotic dependence Vmax ~ Ha is o serve 
including the misaligned one, starting from about B = 0.25 Tesla . 



• i • r- 1 1 r n — Li — i = i 25 • 1 0^ : stream function (a), stream 

Ficrnre 7* Thermal convection without magnetic field for D — 2, li — 1, a« a oc / x 

function for misaligned configuration (b), vertical velocity profils by SEM and FEMR, Vy(x) for (a) at y - 0.25 (c), 
for (b) at v = 0.5 (d), velocity scale is 0.225 cm/s. 


11 

American Institute of Aeronautics and Astronautics 







. d , A n i H 1 Rn - 1 25 • 10 s Ha = 2170(B = 5.0 Tesla)-, stream 

Figure 8: Thermal convection with magnetic field. D - 2, H - l, Ka- k . . velocily 

function (a), same for misaligned case (b), vertical velocity profile Vy(x) for (a) at y = 0.25 (c) and vertical y 

profile Vy(x) for (b) at y = 0.5 (d), velocity scale is 0.225 cm/s. 



Figure 9: Summary of magnetic field suppresion of the flow for H = I. D = 2: maximum value of 
velocity versus S for aligned and misaligned configurations. Predicted asymptotic dependence ^ 

for both cases. 
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THERMAL CONVECTION 


THERMAL CONVECTION 



(a) 


nCRMN. CONVECTION 




(d) 


Figure 10: The velocity profile in the boundary layer (X coordinate from -1.0 to -0 95 ) with ‘"creasing the magne 
field for geometry 2, D = 2, H = 1 , Ra = 1 .25 ■ 1 0 5 : (a) to (d) correspond, respectively toB - 0, 0.5, 5.0 and 5CK0 Tesla. 
The velocity amplitude decreases as Ha and its gradient on the wall decreases as Ha~ with increasing the Hartmann 

number Ha. 
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Abstract 

The Multiquadric Radial Basis Function (MQ) Method is a recent meshless collocation method with 
global basis functions. It was introduced for discretizing partial differential equations (PDEs) by Kansa 
in early nineties. The MQ method was originally used for interpolation of scattered data, and it was 
shown to have exponential convergence for interpolation problems. 

In [11] we have extended the Kansa-MQ method to numerical solution and detection of bifurcations 
in ID and 2D parametrized nonlinear elliptic PDEs. We have found there that the modest size nonlinear 
svstems resulting from the MQ discretization can be efficiently continued by a standard continuation 
software, such as AUTO. We have observed high accuracy with small number of unknowns, as compared 
with most known results from the literature. 

In this paper we formulate an improved Kansa-MQ method with PDE collocation on the boundary 
(MQ PDECB): we add an additional set of nodes (which can lie inside or outside of the domain) a ja- 
cent to the boundary and, correspondingly, add an additional set of collocation equations obtained via 
collocation of the PDE on the boundary. Numerical results axe given that show a considerable improve- 
ment in accuracy of the MQ PDECB method over the Kansa-MQ method, with both methods having 

exponential convergence with essentially the same rates. 

Keywords; Radial basis functions, multiquadric method, numerical solution, continuation, bifurca- 
tions, nonlinear elliptic PDEs. 


1 Introduction. 

The Multiquadric Radial Basis Function (MQ RBF or, simply, MQ) method is a recent meshless collocation 
method with global basis functions, for discretizing PDEs. It was originally proposed in 1970 [19], 120] tor 
interpolation of scattered data and was shown [27], [28], [31] to have an exponential convergence for function 
approximation. The MQ method was introduced for solving PDEs in Kansa [24], [25 m e^ly mneties_ Si c 
then it was successfully applied for solving a number of 2D and 3D PDEs, see e.g. [4], [18], [ J, [ ]» [ by] 
and references there, while some convergence results for solving PDEs, based directly on the interpolation 
error estimates, appeared only recently [14. 15]. Application of the MQ method to PDEs leads to finite 
dimensional problems with full matrices. The Kansa-MQ method was shown to give high accuracy with a 
relativelv small number of unknowns (tens or hundreds for 2D problems). The corresponding linear sys ems 
can be efficiently solved bv direct methods. In [11] we have extended the Kansa-MQ method to numerical 
solution of parametrized nonlinear elliptic PDEs. We presented there results of our numerical experiments 
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with continuation of solutions to and detection of bifurcations in ID and 2D nonlinear elliptic PDEs. We 
found that the modest size nonlinear systems resulting from the MQ discretization can be efficiently contmued 

bv a standard continuation software, such as auto [5]. , 

" Our observations have shown that the residual error is typically largest near the boundary (by one to 

two orders) compared to the residual error in the domain far from the boundary. 

In this paper we formulate an improved Kansa-MQ method with PDE collocation on the boundary 
(PDECB): we add an additional set of nodes (which can lie inside or outside of the domain) adjacent 
to the boundary and, correspondingly, add an additional set of collocation equations obtained via collocation 
of the PDE on the boundary. The motivation for this modification of the Kansa-MQ method comes from 
our observations that 1) the residual is typically the largest near the boundary (by one to two orders larger 
than in the domain far awav from the boundary), and 2) the residual is dramatically reduced when we use 
the PDE collocation on the boundary. The MQ PDECB method leads not only to a higher accuracy, but. 
for nonlinear problems, also to a higher efficiency due to the reduction of the number of unknowns in the 

continuation process by using a preprocessing. . „ 

We apply our MQ PDECB method to several model ID and 2D linear and nonlinear elliptic PDEs 
and present 'results of our numerical experiments. These results demonstrate considerable improvement in 
convergence of the MQ PDECB method over the Kansa-MQ method, with both methods having exponen- 
tial convergence with essentially the same rates. To our knowledge, this is the first demonstration of the 

exponential convergence for the MQ method applied to PDEs. 

A related idea was successfully used for high Re number fluid flows in the cases of the RNS model [8j, 
[12] and Alexeev hydrodynamics equations [10] (in the framework of the finite element method), that was 

applied for the solution of 3D thermo- vibrational flows [9]. . 

\ class of global numerical methods for ID and 2D problems, the numerical algorithms without saturation, 
was proposed by Babenko in early eighties [1]. These include a highly accurate discretization method for 
PDEs based on Chebyshev polynomials. This method was further developed by Belykh (see e.g. [ \,[o\) : 
who found it to be more accurate and better conditioned than the spectral method. 

In Section 2 we formulate the Kansa-MQ and the MQ PDECB methods for a linear elliptic PDE. n 
Section 3 we describe in detail the Kansa-MQ and the MQ PDECB methods for continuation of solutions 
to parametrized nonlinear elliptic PDEs. For clarity of presentation, Section 3 is written independent y o 
Section 2. In Section 4 numerical examples are given that illustrate the accuracy of our method. In Section 
5 we summarize our results. 


2 A linear elliptic PDE. 

We consider a well-posed elliptic boundary value problem: for given g(x), f(x) find u(x) from 

Iu(x) =/(*), in ncW 1 , (1) 

Bu(x)\qq = s( x )i 

where D is a bounded domain. L is a linear elliptic partial differential operator, and D is a boundary operator. 


2.1 The Kansa-MQ method. 


Introduce a set 0/, of nodes (Fig. 1) 

0h = {{*,};=! c n, 


{*i} 


N+N b 
»=N+1 



( 2 ) 


and the MQ basis functions, 

g j {x)=g j (c j ,x) = JWx-x.jWlt+tPj, 3 = h N + N b> g N+ N k+ i(x) = 1, (3) 

where ||x - x.j|| Rd is the Euclidean norm in 3 d . Cj > 0 are called shape parameters [25]. We look for the 
approximate solution Uh to (1) in the form 


;V+N b +l 


(4) 


2 



N b 1 2 . . . N N + 1 

<• © © © •“ 

0 ~ h ~~ 1 


N + N b + 1 N b 1 2 . . . N N + 1 N + 2N b 

: $ (D - © © @ i— 


(a) ID nodes 




+ + + + 


(b) 2D nodes 


+ 

+ 

+ 


+ 


Figure 1: Nodes for the Kansa-MQ and MQ-PDECB methods: (a) ID Kansa-MQ nodes (top) arid PDECB 
(bottom), node numbering is shown; (b) 2D Kansa-MQ nodes (left) and PDECB (right), o no es or 
collocation. . - BC collocation. © - PDE and BC collocation, + - nodes added for PDECB, hi is a distance 
to the boundary (may be negative, if nodes are inside); h is a mean distance between nodes. 


Substituting u h (x) into (1) and using collocation at the nodes 0/,, we obtain the finite dimensional problem 


/.V + iVt+l 


iV+A',, + 1 


L y 52 ajQjixi) j= 52 a J Lg j (x i ) = fi x i), i = 

f,\+,V h +l \ .V+Ni+l 

B ( 52 = E a i B 9 j( x i) = 9 (xi), i = N + + Nb, (°) 

•V+.V fc 

52 aj = o. 

Introducing the notation: a = (<u, ...,rt.\+,v b +i) T ! & = [f{xi),...,f{xff),g{xK+i),...,g(x!s l + l \' b ),0) ^ ® 

Lgi(xi) ... LgN+N b +i( x i) 


L g — 




Bgi(x.\+\) 


Lgi(xx) ... LgK+N b +i{xN) 
Bg l \ + s b {xN+ 1 ) BgN+ff b + i(*jv+i ) 


Boi(iv4-yJ Bgx+.s b {xN+N b ) Bg N +N b +i{xN+N b ) 

l' ... 1 0 


(6) 


, W = 


we can rewrite the system (o) in the matrix form as 

Wa = b, 

whose solution is 

a = 1 V~ l b. 


(7) 

(8) 


2.2 The MQ PDECB method. 

Introduce a set 0/, of nodes (see Fig. 1) 


0, 


= {{x,}; N =l c n, {x.iEv+i c an, {x t }f = + 2 ;\ +1 c \ an} 


(9) 


3 




where the nodes {*,},=^+* l+1 , which can be inside ft or outside ft, axe adjacent to the boundary 9ft, and 
the MQ basis functions, 


gj(x) = Oj(cj,x) = y/\]x - ZjIIr, + Cj, j = l,-,N + N b ,g N+ 2 N t+ i(x)-l. (10) 

We look for the approximate solution u h to (1) in the form 


N+2iV 6 + l 

Uh(x) — E a,j9j{x), 


( 11 ) 


Substituting u^(x) into (1) and using collocation at the nodes 0^, we obtain the finite dimensional problem 


( A r +2A r b+l \ 

E aj<7 ; (xt) j = E djLgj(xi) = f{xi), z = 1, A r -h Ae>, 

/A'+2A r b + l \ A r +2A r h + l 

2? = E i = iV -h 1, - A r + A&, 


/V+2A r b +l 


( 12 ) 


j=i 


j=i 


A>2A r h 


Y a i =0 - 


j=l 

Introducing the notation: a = (ai , aAr+2A : b +i) T ? ft = (/( x i)> /(^iV+A^i)* 5( x N+iVb+i)i 0( x ;V+2A b ), 0) € 

*^A* + 2 A\ -j-1 


Lg — 


Lgi(xi) 


FffA , +2,v b +i(£i) 


Z? 3 = 


L Lgi(x\+\ b ) I-3/V+2,Vi,+i(a;jV+N b ) 

Br?i(x,v+i) ••• Bg N +2S\{xN+i) Bg N +2N h +i{xN+i) 


Bfh{x.\+ \ t 

1 


£?g,v+2.V h (x,V + ,V b ) S5JV+2,Y b + l( I N+A' b ) 
1 0 


(13) 


, W = 




we can rewrite the system (12) in the matrix form as 

XV n — /l 


d at 


3 Continuation for nonlinear elliptic PDEs. 

Consider a boundary value problem for a second order system of n parametrized nonlinear elliptic PDEs: 


F (u(x). A) = B(a)Au - /(Vu, u, x, A) = 0, in ft C K d , A € K, ti(-) € S n • (15) 

2 ?u(x)|^q = 0 . 


where ft is a bounded domain. D{ X) is a positive diagonal n x n matrix, / is smooth, and i? is a boundary 
operator which we assume, for simplicity, to be linear. For the bifurcation analysis m the process s of con- 
tinuation we also need to consider the eigenvalue problem for the linearization D iF(u,A) of F about the 

solution u of (15) 


D { F[u , \)v(x) = /iv(x), in (1.6) 
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3.1 The Kansa-MQ method. 

To formulate the approximate problem, we first introduce the set 0/t of nodes 

e, = {(*,}", c n, {*,}£& can} 

and the MQ basis functions, 


9j (x) =gi(Cj,x) = + j = \,...,N + N b , g N+Nb+l {x) - l, 

We next define an MQ finite dimensional subspace 

{ M+N b +l s+N, . • . • „ .. 

X= Z »<»(•) > E “1 = °' b>k*.) -o,t-n + v .., N + i», 

;=1 1=1 

The problems (15) and (16), respectively, are approximated by the collocation equations 

F (uh(x,), A) = 0, u h €S h ,i = 1,...,1V, 

Lt» h (*i) = D t F(u h , A)v h (Xi) = ttv h (Xi), v h € S fc , i = 


Substituting 


iV + Nb + l 

Uh(x) = Y a 3»j( x )- 

J =1 

iV + ^b + l 

t’h(x) = Y W x )> 

j=i 


(17) 

(18) 

(19) 

(20) 
( 21 ) 

( 22 ) 

(23) 


into (20) and (21), respectively, and using the definition (19) of S k , we obtain the following finite dimensional 
problems: 

/n-hv b +i \ 

(G(o.A)),=fI Y a }0j(xi)’ X J =0 ’ t = 1, JV, 


//V'+,Y b -rl \ 

B ( Y a i 9 iMj = 0 ’ i = N + l,...,N + N b , 

N+Nb 

Y a > = °* 

j=i 


(24) 


( iV+iVb + l 

( N+Nb+1 
jV+ A r b 

E ^ = 0 - 

j=i 


JV+iV b + l 

= n E ^(*0. z = i, JV, 

j= 1 

= 0, t = A r + l,...,iV + N b , 


(25) 
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X T i /i k \T fz TIJ n x ( TV + N’t, + 1 ) 

Introducing the notation: o = (ai, ...,aN+\ b -ri) > & — v°i > •••’ °N+N b +i) 

J3<7i(z/v+i) ••• B<?;v+w b (x,v+i) Bgw+w b +i( x /v+i) 

9 Bgi{xy+N t ) ■■■ Bg^+N b { x N+N b ) BgN+N b +l{ x N+N b ) 

1 ... 1 o J 

fll(zi) — 9N+N b + l{ x l) 
gi{XN) 9N+N b +l( x N ) J 

we can rewrite the problems (24) and (25) in the matrix form as 

G(a, A) = 0, 

B g a = 0, 

L g b = nTb, 

B g b = 0. 

Implementation 1. Let 

a 1 = (ai,...,u/v) T € R N , a' = {a N+ u -,aN+N b + 1) G K 



' Lgi(xi) ... 

Lgs(x i) 


Lff.v+i( x i) 


L\ = 




- = 


• 

9 

L L 9i {xn) - 

Lg N {xs 


L<7,v+i(xn) 

... LgN+N b+ i(Xjv) 



' Si( x i) 

... 3,v(x 1 


’ 9.v+i( x i) -- 

• SfA'+jVfc+l (*^ L ) 

r 1 


| !7i( x iv) 

... 9n{x.\ 

, r 2 = 

. 3.v+i(xat) •• 

■ 5A'+iV b +l(*iv) . 


Bc\ (x/v+i) - Bg\(x N+ 1) 

Be: (iN+ivJ - Bg.v(x N +N b ) 

1 ... 1 




L. 

F Lffi( x i) 

... Lg\ + .\\^\ (x, ) 


^ 9 = 



, r = 


Loi (x,v) 

... LgN+s\-ri (Xjv) . 



Bp.v+i(xn+i) ••• B(7/v+A' b (x.v+i) Bgjv+iV b +i( x /v+i) 

3 B3v+i(XjV+N b ) ••• B5jv+,v b (x,v+iV b ) BgN +1 \ b +i{xN+N b ) 

_ 1 ... 1 0 

Substituting this into (27), we rewrite it as: 

g:a\\) =G{a l ,a 2 .X) =0, 

where a 2 solves 


( 26 ) 


(27) 


(28) 


(29) 


(30) 
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Similarly, we rewrite (28) as 


L l g b l + L 2 g b 2 = n{T'b l +T 2 b 2 ) , 
B l g b l + B^b 2 = 0. 


or, eliminating fr, as 

L l g b l - L] (B 2 )" 1 Bib 1 = M (rV + r 2 (B 2 )' 1 B]b l ) . 

We are interested in continuation of solutions to (29). Therefore, in addition to a 1 
unknown, and add an algebraic constraint 

GcCo 1 , A) — 0,' 


(31) 

we also treat A as 


(32) 


which defines a parametrization of the solution curve. . . 

Algorithm 1 (Continuation algorithm for the system (29), (32)). Given curren approxim i 
a 1 € S nx ' ,v and UK, a complete Newton iteration consists of the following steps: 


( 0 ) Compute the matrices B’, B 2 , F 1 , T 2 . 

( 1 ) Solve the system (30) to find a 2 . 

( 2 ) Use the expressions (29), (32) to compute the residuals -G(a\ A), -G c (o\A) and then compute the 

matrices D iQ, D2Q1 D\G C , and D2G C by differencing. 

(3) Solve the system 

D X C8 a 1 +8\ D 2 Q = -G{a\ A), ( 33 ) 

DiG^a 1 + 8XD2GC = — G c (u x ,A), 


where we omitted iteration indices for 5a l and 5 A in (33). 

( 4 ) Update a 1 — ► a 1 4* tia 1 and A -4 A 4- JA. 

(5) Solve the generalized eigenvalue problem 

DiQb 1 = y. (r 1 + T 2 (B 2 ) -1 Bg) b l (3-i) 

(to detect bifurcations). Note that D\Qb 1 = L l g b — L g (B g ) B g b , see (31). 

Implementation 2. Let U = (V\, Us T be the vector of nodal values of the solution u„ (22) of the 
collocation problem (20), and let (<M£i> be the Lagrange basis m S h : 

{^€S h :Oj(i,)=^ = (35) 


Then Uh can be written as 


N 


Uh(l) — 


}= 1 


Combining this with the definitions (22) of Uh and (26) of B g and T, we have 


ra = U, 
B g a = 0, 


which defines the 1 - 1 correspondence between U € 


and a € !R"x(A'+v k +i). 


(36) 


(37) 
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( 38 ) 


The problems (20) and (21), respectively, are written as 

g(U,X) =G{a,X) =0, 


where a solves (37), and 

D X Q{U.\)V = txV, V € R nxiV . 

As before, to define a parametrization of the solution curve, we add an algebraic constraint 


G c {U,\) = 0. 


(39) 

(40) 


Algorithm 2 (Continuation algorithm for the system (38), (40)). Given current approximations to 
U e K nx ' v and A e K, a complete Newton iteration consists of the following steps. 


(0) Compute the matrices B g ,T. 

(1) Solve the system (37) to find a. 

encing. 

(3) Solve the system 


D x Q5U + 6XD 2 G = -G(U, A), 
DxG c 5U + SXDiGc = - Gc{U,X), 


(41) 


where we omitted iteration indices in (41) for SU and <5A. 

(4) Update U -* U + SU and A -* A + <5A. 

(5) Solve the eigenvalue problem (39) (to detect bifurcations). 

Remark 1 For our numerical experiments, we implemented in AUTO [ 5 ] Algorithm 2 for sa-MQ 

no‘ ,S Zt e Zr J* \L **«- f W > . 

generalized eigenvalue problem whose solution is not supported by AUTO. 


(42) 


3.2 The MQ PDECB method. 

To formulate the approximate problem, we first introduce the set Q h of nodes 

where the nodes {x,)£S i,. t , , which can be inside a or outside H, are adjacent to the boundary Ml, and 
the MQ basis functions, 

gj(x) S gj(Cj,x) = yj\\x - *,;& + <?. i = l,-,^ + dY 6 J ff/v+ 2 A l+ i(x) = l. (43) 

We next define an MQ finite dimensional set (which is not a subspace, in general) 


{ jV+2.V b + l 

x— E a ’ 9 ^ ] : 


N+ 2.V b 

E 

j=i 


2 v b 

aj =0, Bx(x0 =0, F(x(x 4 ),A) =0, i = N + l,...,iV + iV b 


(44) 
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The problems (15) and (16), respectively, are approximated by the collocation equations 

F (uh(xi), A) =0, u h e Si i = 1, JV, 


Lvh(Xi) = DiF(iifc,A)fh(Zi)=M«fc(*i)i = 


Substituting 


iV+27V b + l 

Uh(x) — 

i =1 

jV+27Vb+l 

Wh(*) = 

3 = 1 


(45) 

(46) 

(47) 

(48) 


into (45) and (46), respectively, and using the definition (44) of S£, we obtain the following finite dimensional 
problems: 

(G(«,A)), = F feS 2iVt+1 =0, i = 1,-.*. 

(G(o,A)) j HF(EfJ; aAr ‘ +1 Ojfli(*0,A) =°, » = jV + l,...,JV + iV 6 , (49 ) 

B SS SA ' +1 «jSj(*0)=O, i = rY + l,...,N + JV b , 


ES aJVl oi=o, 


' JV+2JW + 1 


jV+2JV b + l 

E 

j=i 


i( E «*<>)=" £ * = l,...,iv + isrfc, 

( i\ r +2.V b + l 

B 


Y bjOjixi) 1=0, * = N + 1 ,...,N + N b , 

V / 

.V-K2A\ 

E b ; =0 - 

i=i 

T i ru k \T ^ ranx(JV+2A r 6 + l) 

Introducing the notation: a = (rti, 0 JV+ 2 /v b +i) » ^ /v>2/v b +u 

r B< 7 i(x. v+ i) ... Bg N+2 s b {x s+ i) Bg N + 2 N b +i{xN+i) 


B 0 = 


L 3 = 


Bgi{xN+Ni) 

1 

£,<7i(xi) 


. . B</w + 2A' l (i.v+.v b ) Bgjv+2jv k +i(a:jv+jv fc ) 


(50) 


(51) 


1 0 

J 


£tfA r +2N b +l( x l) 

, r = 

■ <7i(zi) 

. 0N+2jV b + l( x l) 

EtfN+2;Y b +i( x .v-i-Ar b ) . 


_ 9 l(XN+N b ) ■■ 

■ 5N+2A r b -rl( x .V+N b ) . 


Li 7i(x.y+,vJ 

we can rewrite the problems (49) and (50) in the matrix form as 

(G(o. A)),=0, t=l,.-.,iV, 

(G(a,A)), =0, i = xV + 1,-,N + 1V 6 , 
B g a = 0, 


L 9 b = /xrb, 
B g b = 0. 


(52) 


( 53 ) 
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Implementation 2a. Let U = (u h (zi), ...,u h (x N )) T be the vector of nodal values of the approximate 
solution u h . Then by the definitions (47) of u h and (51) of B g and V, we have 


(G(a,\))t = 0, i = N + \,...,N + N b , 

Ta = U, 

Bgll = 0, 

which defines the 1 - 1 correspondence between U eW'* N and a e R nx{N+ 

The problems (45) and (46), respectively, are written as 

(Q(U t \))i = {G(a, \))i = 0, i = h-,N, 

where a solves (54), and , * 

DiG(U, A)V = mV VeR n * N . 

As before, to define a parametrization of the solution curve, we add an algebraic constraint 

G C [U, A) = 0. 


(54) 


oo 


(56) 


(57) 


Algorithm 2a (Continuation algorithm for the system (55), (57)). Given current approximations to 
U e RnxJV and A e 3, a complete Newton iteration consists of the following steps. 

(0) Compute the matrices B g , T. 

(1) Solve the system (54) to find a. 

(2) Use the expressions (55), (57) to compute the residuals -SffJ , ,»), “n C^TlTy differ- 

matrices D,f • D.S(C.A), D,9 m D.& = M - D,9AV,» b> **» 

encing. 


(3) Solve the system 


DiGSU + 6\D 2 g = -G{U,\), 
DiG c fiU + SXD 2 G C — ~ Gc[U,X), 


(58) 


where we omitted iteration indices in (58) for 6U and 6\. 

(4) Update V -t l T + 6U and A -t A + SX. 

(5) Solve the eigenvalue problem (56) (to detect bifurcations). 

4 Numerical experiments for ID and 2D elliptic PDEs. 

We present examples of solution of linear ID and 2D elliptic PDEs and continuation of solutions to nonlinear 
m and 2D GeSd-Bratu equation. Each problem is discretized by the Kansa-MQ method, see Eq. (38), 

and tatfit! cilof^onlinear pmbTem^wepwform continuation of solutions by Algorithm 2 ^theKansa^MQ 
method and bv Algorithm 2a for the MQ PDECB method. We compare the accuracy of the detection 
the limit point (or fold) bv the two methods. We recall that a solution (u 0 ,A 0 ) of equation /(u A - 0 
U a SS limit point if 'the solution curve in („(,) A(s)) for some par^etrization s, makes a turn at 

(uq. Aq). This is expressed formally as dim A {f u { u o,^o)) 811 A( u 0i o) ' $ v tx °_L , _ n(is _ \\ 

We will use throughout the notation h for the average distance between the nodes. Then h - l/[K ) 
for a ID problm on (0,1) and for a 2D problem on (0,1) x (0,1), where K is the number of nodes along 

eaL To improve the accuracy, we employ 2 simple adaptation strategies for the shape ^ “ 

{C1 c .v+.v b } for the Kansa-MQ method, see Eq. (18), and C 1 = (c 1 , c.v +2 A h yor the MQ PDECB 

method.' see Eq. (43); for the nodes Q h for the Kansa-MQ method, see Eq. (17), and 0, for the MQ PDEC 
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method see Eq. (42). To be specific, assume that Q = (0, 1) x (0, 1) and consider the case of the Kansa MQ 
method’ Let r(x,y,C,0 h ) be the residual. Our strategies are all based on the Nonlinear Least Squares 
Method which minimizes the L 2 norm y>(C, Q h ) = ||r|| 2 of the residual. By the quasx-umform distribution 
of nodes we will mean the distribution of nodes, where the nodes adjacent to the boundary W are placed 
at the distance h = 5h 0 , Q<8<\, from 30. while the remaining nodes are distributed uniformly with the 

distance h 0 between them. 


Strategy 1 . Uniform distribution of nodes 0/,; Ci — ... — — c; miny?(C, Oh), 

Strategy 2. Quasiuniform distribution of nodes 0j,; c\ = ... = CN+N b = c \ min y>(C, 0h). 


In all examples below we use the adaptation strategy 2.. . - 

Example 2 .4 ID model linear problem 

u xx + (2")' sin(27rx) =0, in fl = (0, 1), 
u(0) = u(l) = 0. 


The analytical solution is 


^ex act — sin(2nx). 


Numerical results are presented in Fig. 2a. 


(59) 


ID PDE BY PDECB AND KANSA-MQ 



(a) ID linear PDE solution 


ID BRATU-GELFAND EQUATION BY PDECB 



(b) ID continuation 


Figure 2: Convergence properties of the Kansa-MQ method and the MQ PDECB method. 

(a) ID linear problem, Eq. (59); the L x norm of the solution error is plotted, in the logarithmic scale versus 

1 jh where h is the average distance between the nodes. The roundoff error starts to domina e a / ~ 

for Kansa-MQ method and at l/h « 18 for the MQ PDECB method. 

(b) The location A of the limit point for ID Bratu-Gelfand problem, Eq. (61). Relative error in A is plotted 
in the logarithmic scale versus l/h. 


Example 3 A 2D model linear problem 

Au-(2xV+2x 2 y + 2xr-6xy)e( I+ '' ) )=0 in fl = (0, 1) x (0, 1), ( 60 ) 

u 0 . 
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(a) 2D linear PDE solution and interpolation 


ID BRATU-GELFANO EQUATION 



Fie-ure 3- Convergence properties of the Kansa-MQ and the MQ PDECB methods. 

fa? 2D linear PDE, Eq. (60); the L TC norm of the solution error is plotted, in the logarithmic scale versus 
.//.where h is .he average distance between the nodes. The roundoff error starts to dommate at 1/1. ~ 9 
for the K^sa-MQ method and at 1/1. ~ 11 for the MQ PDECB method. We also prov.de, for companson, 

thp prror in the MQ interpolation of the exact solution u exa ct- . , 

The location A of the Lit point for 2D Bratu-Gelfand problem, Eq. 62). Relat.ve error m A .s plotted. 

in the logarithmic scale, versus 1 /h. 


The analytical solution is 

U' Xa: ;=x(x-l)y(y-l)e {z+y) . 


Numerical results are presented in Fig. 3a. We do not have an explanation of why the MQ PDECB solution 
is more accurate them the interpolation. 


Example 4 ID Gelfand-Bratu problem. This is a scalar problem 


u" — Xe u = 0, in fi = (0,1), 
u(0) = u(l) = 0. 


(61) 


that appears in combustion theory and is used as the demo example exp in AUT097 [5](fifthorder adaptive 

~ < , e ** 

point is shown in Fig. 2b. See also [11] for additional numerical results and references. 


Example 5 2D Gelfand-Bratu problem 


Au + \e J =0, in n = (0,1) x (0,1), 
u |ar.= 0. 


(62) 


Tins problem was studied by a number of authors. In [6] the problem 

ortho, onal spline collocation method with sparse Jacob, an There is a limit (fold) pc ml tm the sahl 
cure, 77.e exact location of the limit point is assumed to be at the value o A ohta.ned n «/ on 16 x 16 
mesh with 4x4 collocation points. The relative error in location of the limit point is shown in F g. . - 

that the curve for the Kansa-MQ method was obtained [11] using also 

slowed down computations, while we use only double precision with the MQ PDECB method here. 

[Ill for additional numerical results, references and a discussion of the operation count. 
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Figure 4: (a) ID linear problem with a boundary layer, Eq. (63) with e 10 ■ The MQ PDECB solutio 

with 21 nodes and the analytical solution u exact are plotted vers^ x in ( , -b PDECB 

(b) The residuals for the solutions of ID Gelfand-Bratu problem (61) by the Kansa-MQ 
methods axe plotted versus x with 1 /h = 9. A = 2.5. The L x residual norms are 4.2 x 10 and 3.3 

respectively. 


Example 6 .4 ID model linear singular perturbation problem studied m [23]: 

etij, 4- u x = 0, in D = (0, 1), 
u(0) = 0, u(l) = 1. 


(63) 


The analytical solution is 

IW = ( l-e- l/e )/(l-e- 1/f ). 

It was demonstrated by Hon that this problem can not be solved by a standard Kansa-MQ approach for e « 1, 
and the adavtive technique proposed can be an efficient way to treat such problems / 23] ■ 

Here wTlZZeMQ PDECB to solve this problem for, « l with rclotoely small number of nodes. For 
, - 10" 3 the MQ PDECB solution with 21 nodes and the exact solution are plotted versus x in ig. 4_ 
The L no™ of the solution error is 0.001 for the MQ PDECB method, while it is 0.22 for the Kansa-MQ 
mZhJwUh 101 nodes (not shown). Note that for e = 10- one con niton the same error 0.001 .n the MQ 
PDECB solution with 41 nodes (not shown). 

Example 7 Fig. 4b shows the residual distribution for the solutions of 1 D Gelfand-Bratu problem (61) by 

TZsa-MQ Id the MQ PDECB methods with 1/h = 9 A = 2.5. 

and 3.3 x 10" 6 , and the L 2 residual norms are 7.8 x 10 and 1.1 x 10 for the Kansa MQ ana 
PDECB methods , respectively. 

5 Conclusions. 

PDECB) Ttekteotl ,f“e (dtoaba ot'notelljLnt to thebou, clary tmd, corre- 

” ilh b » th methods havine 

exponential convergence with essentially the same rates. 


13 



Acknowledgments 

This work was supported in part by National Aeronautics and Space Administration through jgr a»t ^NAGS- 
1229. The first author is grateful to V.N. Belykh for fruitful discussions, that led us to the idea of the 

PDECB. 


References 

[1] K.I. Babenko, Foundations of numerical analysis, Moscow, Nauka, 1986, 743 pp. (in Russian). 

[2] V.N. Belykh, Algorithms without saturation in an axially symmetric boundary value problems, Soviet 
Math. Dokl . , 36, no.l (1988), 146-150. 

[31 V.N. Belykh, Overconverging unsaturated algorithms of numerical solving the Laplace equation // Math- 
ematics in Applications. Proc. of International Conference honoring academician Sergei K. Godunov. 
Novosibirsk, 1999, (to be published). 

r 4 l G B Cook, M.W. Choptnik. M.R. Dubai. S.Klasky, R.A. Matzner, and S.R. Oliveira Three dimensional 
initial value data for collisions of two black holes, in Approaches to Numerical Relativity, editor R. 
d’Inverno, Cambridge University Press (1993), 265-285. 

[51 E.J. Doedel, A.R. Champneys, T.F. Fairgrieve, Yu.A. Kuznetsov, B. Sandstede, ^ nd X ^' 

AUT097: Continuation and bifurcation software for ordinary differential equa ions (wi om 

(1997). 

[61 E J Doedel H. Sharifi, Collocation Methods for Continuation Problems in Nonlinear Elliptic PDEs, 
Proc. ERCOFTAC and EUROMECH Colloquium, Continuation Methods m Fluid Mechanics, Aussois, 
France, September 1998, to appear in Notes on Numer. Fluid. Mech. 

{71 G.E. Fasshauer, Solving differential equations with radial basis functions: multilevel methods and 
smoothing, Adv. Comp . Math. 11, No. 2/3 (1999), 139-159. 

[8] A. I. Fedosevev, A regularization approach to solving boundary layer problems for Navier-Stokes equa- 
tions. CFD Journal, 11, No. 1 (2000), (in print). 

[91 A. I. Fedosevev and J.I.D. Alexander, Investigation of Vibrational Control of Convective Flows in 
Bridgman Melt Growth Configurations. J. of Crystal Growth, 2000 (to appear). 

[101 A I Fedosevev B. V. Alexeev. Higher Order Continuum Model for Incompressible Viscous Fluid Flow. 
11 “iwZ ltaW. .hi Discrete Systems, Eds. E. Inha, A. Markov, World Sclent, fic, London, 1998, 

130-137. 

[11] A. I. Fedosevev, M.J. Friedman, and E.J. Kansa, Continuation for 

Equations Discretized by the Multiquadric Method, Preprint ^ a ^ A /9812013 a ^ E-PMNT LAN 
(http://xxx.lanl.gov/ps/math.NA/9812013); to appear in Int. J. Bifur. & Chaos 10, No. (2000). 

[121 4. 1 Fedosevev. E.J. Kansa, C. Marin, and A.G. Ostrogorsky, Magnetic field suppression of senucon- 
J ductor melt flow in crystal growth: comparison of three methods for numerical modeling, CFD Journal, 

11. No. 1 (2000), (in print). 

[13] R. Franke, Scattered data interpolation: tests of some methods, Math.Comp., 38 (1982), 181-199. 

[14] C. Franke and R. Schaback, Solving partial differential equations by collocation using radial basis 
functions, Appl Math. Comp., 93, No. 1 (1998), 73-91. 

[15] C. Franke and R. Schaback, Convergence order estimates of meshless collocation methods using radial 
basis functions (1998), preprint. 

[16] M.A. Golberg, and C.S. Chen, The theory of radial basis functions used in the dual reciprocity boundary 
element method, Appl. Math. Comp.. 60 (1994) 125-136. 


14 


[17] M.A. Golberg and C.S. Chen, Improved multiquadric approximation for partial differential equations, 
Engin. Anal with Bound. Elem. 18 (1996), 9-17. 

[18] M.A. Golberg and C.S. Chen, Discrete projection methods for integral equations , Comput. Mech. Publ., 
Boston, MA (1997). 

[19] R.L. Hardy, Multiquadric equations of topography and other irregular surfaces, J. Geophys. Res., 76 
(1971), 1905-1915. 

[20] R L. Hardv, Theorv and applications of the multiquadric- biharmonic method: 20 years of discovery, 
Comput. Math. Applic., 19, No. 8/9 (1990), 163-208. 


[21] Y.C. Hon and X.Z. Mao, An efficient numerical scheme for Burgers equation, Appl.Math. Comp., 95 

(1998), 37-50. . ' • • 

[22] Y.C. Hon, and R. Schaback, On ansymmetric collocation equations by Radial Basis Functions, preprint 
(1999). 

[23] Y.C. Hon, Multiquadric collocation method with adaptive technique for problems with boundary layer, 
Int. J. Appl. Sci. Comput., 6, No. 3 (1999). 


[24] E.J. Kansa, Multiquadrics-a scattered data approximation scheme with applications to computational 
fluid dynamics-I. Surface approximations and partial derivative estimates, Comput. Math. Applic., 

No. 8/9 (1990), 127-145. 


[25] 


E J Kansa, Multiquadrics-a scattered data approximation scheme with applications to computation 
fluid dynamics-II. Solutions to hyperbolic, parabolic, and elliptic partial differential equations, Comput. 
Math. Applic., 19, No.8/9 (1990), 147-161. 


[26] E.J. Kansa, Y.C. Hon, Circumventing the ill-conditioning problem with Multiquadrics radial basis 
functions: Applications to elliptic partial differential equations, Adv. Comp. Math. (1998). 


[27] W.R. Madych, Miscellaneous error bounds for multiquadric and related interpolants, Comput. Math. 
Applic. 24, No. 12 (1992), 121-138. 

[28] W.R. Madych and S.A. Nelson, Multivariate interpolation and conditionally positive definite functions 
II. Math. Comp., 54 (1990), 211-230. 


[29] M. Sharan, E.J. Kansa, and S. Gupta, Applications of the multiquadric method for the solution of 
elliptic partial differential equations, Appl. Math. & Comput., 84 (1997), 27o-302. 

[30] Wong, S.M., Y.C. Hon, T.S. Li, S.L. Chung, and E.J. Kansa, Multi-zone decomposition of time- 
dependent problems using the multiquadric scheme, to appear m Math. Comp. Applic. (1998). 


[31] Z. Wu and R. Shaback, Local error estimates for radial basis function interpolation of scattered data 
IMA J. Num. Anal. 13 (1993), 13-27. 


15 


Computational Fluid Dynamics Journal, 
vo 1 . 9 No . 1 (2000), (in print) 


A REGULARIZATION APPROACH TO SOLVING THE NAVIER-STOKES 
EQUATIONS FOR PROBLEMS WITH BOUNDARY LAYER 


Alexandre I. Fedoseyev 


Center for Microgravity and Materials Research, 
University of Alabama in Huntsville, 
Huntsville, Alabama 35899, USA 
E-mail: alex@cmmr.uah.edu 
http://uahtitan.uah.edu/alex/ 


Abstract 

We present an alternative method for treating numeri- 
cally a problem of viscous flow with a boundary layer 
that is based upon regularization of the Navier-Stokes 
equations. Shishkin (1997) showed that grid methods 
perform poorly in dealing with the boundary layer. Tra- 
ditional grid methods give poor agreement with the ex- 
perimental data for high Re number flows. Shishkin 
showed that the remedy for these difficulties is the con- 
struction of special meshes in boundary layer. We 
present an alternative approach that is more efficient 
and less mesh dependent. Our approach is based upon 
the regularization of the Navier-Stokes equations, and 
we discuss the mathematical and physical aspects of 
this approach. Numerical results that we obtained by 
our regularization process in 2D and 3D are compared 
with the experimental measurements. We compared our 
model against : (1) The 3D driven cavity flow by Kos- 
eff and Street (1982) at Re = 3200 and 10.000; (2) a 
2D backward facing step flow by Kim et al. (1980) 
at Re = 44.000; and (3) a 3D thermal convection in a 
cylinder by Bogatirev et al. (1996) at Ra = 1000 to 
100.000. This proposed regularization model is not a 
turbulence model, and no additional equations are in- 
troduced. Recipes for the choice of the regularization 
parameter are presented. 

Keywords; Navier - Stokes equations, Boundary Layer, 
Regularization, Finite element method 

1 Introduction 

Driven cavity problem is widely used as a benchmark 
for comparison of numerical codes. The agreement be- 
tween different codes is within 1% or better. Published 
2D Navier-Stokes (NS) solutions can qualitatively de- 


scribe the flow structure, the number and location of 
vortexes and their size, but show poor agreement with 
the experimentally measured velocity profiles by Kos- 
eff and Street (1984) for Re = 3200 and 10 4 (see [1] 
and references therein). 

Results obtained by Ghia et al. [2] on a fine mesh 
(256x256) for Reynolds number up to Re = 10 4 are only 
the stationary solutions. While the actual fluid flow is 
essentially transient and 3D, the measured mean veloc- 
ities in the plane of symmetry (y = 1.5) appear to be 
2D. A disagreement of 2D numerical solution with ex- 
perimental data is by a factor of 2 to 3 (Fig. 1, 2). A 
3D NS solution still can not improve this discrepancy, 
and known 3D results differ significantly (see 1992 
GAMM-workshop [3]). One may justify that the flow 
turbulence is a reason for a disagreement. We tested 
this and solved the problem using a standard k-E turbu- 
lence model that is available in most commercial CFD 
codes. Our results obtained with the CFD2000 code [4] 
are presented in Fig. 1, 2 as well. One can see that the 
agreement of the k — £ model results with experiments 
is still very poor. 

We observed a similar discrepancy with the exper- 
iment in the numerical modeling of weakly turbulent 
thermal convection in a vertical slot (H/L = 11.2) 
heated from a side, the Rayleigh number was Ra L = 
3.75- 10 8 , and Pr = 15. A direct numerical solution (us- 
ing Boussinesq approximation) allowed us to describe 
basic features of the flow evolution and large vortex dy- 
namics [5],[6]. 

Mean vertical temperature profile obtained agrees 
well with the experimental data [7], but mean velocity 
profiles differ significantly from those data. Local mag- 
nitudes of computed velocity were about twice higher 
than ones in the experiments [7]. The Reynolds number 
based on the computed velocity was of the order 10 4 . A 
use of more detailed time and space discretization did 
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not improve the results. 


Re=3200 HORIZONTAL VELOCITY PROFILES 



Figure l: Driven cavity problem, Re = 3200: Compar- 
ison of horizontal velocity profiles at a = 0.5 for nu- 
merical (lines) and experimental (squares) results . 1- 
Navier-Stokes solution, 2 - k — £ model solution 

Recently a poor ability of grid methods to deal 
with boundary layer has been proved theoretically by 
Shishkin [S]. The estimation for the solution error is 
given as 0(1) for uniform meshes, if thin boundary 
layer is present. Shishkin pointed out that the rem- 
edy for these difficulties is the construction of a special 
meshes in the boundary layer (now called as Shishkin 
meshes). To construct a mesh, one needs to use the 
boundary layer thickness, that is not known in advance. 

Shishkin's theory received confirmation for a driven 
cavitv llow problem in [9]. It was shown that the 
computational results are extremely mesh sensitive for 
Re - 3200 and higher. Still, with some special fine 
mesh resolution near the walls, proposed in this paper, 
it was possible to obtain a solution that agrees well with 
the experimental data [1]. 

In this paper, we present an alternative approach that 
is efficient and less mesh dependent. It is based on the 
regularization of the Navier-Stokes equations. Mathe- 
matical and physical aspects of this approach are dis- 
cussed. Presented model is a major simplification of the 
higher order continuum model for incompressible vis- 
cous flow proposed in [ 10]. 

The governing equations are presented, a model is de- 
scribed' and applied to a several viscous fluid flow and 
thermal convection problems. Comparison with the ex- 
perimental data and results obtained by different meth- 


Re«3200 VERTICAL VELOCITY PROFILE 



Figure 2: Driven cavity problem, Re = 3200: compari- 
son of vertical velocity profiles at z = 0.5 for numerical 
(lines) and experimental (squares) results : l -Navier- 
Stokes solution, 2 - k — £ model solution 

ods are discussed as well. 

2 Governing equations 

We consider a flow of incompressible viscous fluid 
in a closed domain. The governing equations are the 
following. The momentum equation is 

^X + (VV)V-/?e“ 1 V 2 V + Vp - F = 0 (1) 

3t 

while continuity equation is 

V ■ V = 0 (2) 

where Re = VqL/v - the Reynolds number, V'o - velocity 
scale, L - hydrodynamic length scale, v - kinematic vis- 
cosity, Fisa body force. For a case of thermal convec- 
tion the body force is F = GrRe~ [ ■ 0 • e K (Boussinesq 
approximation), where 0 is a nondimensional tempera- 
ture, Gr is the Grashoff number and is a unit vector 
in the direction of gravity. In the case of thermal con- 
vection the energy equation complete the formulation: 

— + (\V)Q = Pr~'Re-'V 2 Q (3) 

where 9 is a nondimensional temperature scaled by 
0 = (7 - T cnl< i)/&T. AT = T hll , - T al ,j . The Prandtl, 
Grashoff and Rayleigh numbers are respectively Pr = 


2 




v / k, Gr = Ra/Pr , and Ra = $ATgL*k~ l v where j3, 
g, k are respectively the coefficient of thermal expan- 
sion, gravitational acceleration, and the thermodittusiv- 

ity. 

Boundary and initial conditions accomplish the prob- 
lem formulation and will be given for each of problem 
considered below. 


The value of T (or /) is not known in advance, and we 
provide recipes on a choice of these values at the dis- 
cussion of the results. 

Such a regularization has an additional useful feature 
for the FEM. It allows to use the same order finite ele- 
ment approximation for the velocity and pressure with 
all the unknowns located at the same nodal points. 


3 Regularized Navier - Stokes 
equations 

The regularization proposed consist in the modification 
of the continuity equation (2) that becomes as follows 

V • V = x“ V • (Vp — F) (4) 

where x* is a small regularization parameter. For x y 
0 we approach the original equation (2). A boundary 
condition for the pressure on a wall is 

(Vp — F) • n = 0, (5) 

where n is a unit wall normal vector. Equations (4) 
and (5) present the main feature of this method, and en- 
sure a balance of the normal component ot the momen- 
tum on the wall. We will show in the following sec- 
tions that equations ( 1 ), (3), and (4) with (5), called here 
as regularized Navier-Stokes equations (RNS), are more 
suitable for use with the finite element method ^FEM). 
The solution of the RNS by FEM can result in a bet- 
ter agreement with the experimental measurements for 
high Reynolds number flows. The RNS model showed 
also a good agreement with the analytical solutions for 
a very slow flows with thin boundary layer in the case 
of magnetic field suppression of the semiconductor melt 
flow in crystal growth [11]. 

For a justification of this regularization one can be 
referenced to the new hydrodynamic equations pro- 
posed in [ 1 2] that have the fluctuation terms of this kind. 
Equation (5), according to the theory in [12], means the 
absence of the hydrodynamic fluctuations on the wall. 

Similar term, as in the right-hand side of eq. 4), has 
been obtained as a result of the “consistent treatment 
of time-advancement for the divergence-equation by 
Ldhner [13]. A discretization time step A/ was used in 
[ 1 3] instead of parameter x ¥ as in our RNS model. 

A nondimensional regularization parameter X* is 
expressed through dimensional scales as following, 
t * — tL^Vo, where x is some time scale. Note, that 


4 Problems solved 

We apply the RNS model for the solution of the follow- 
ing problems: 

1) Lid-driven cavity problem (2D and 3D). We 
compare our results with experimental data obtained by 
Koseff and Street[l] for Re - 3200 and 10 4 . We solve 
equations (1) and (4) numerically for the unknowns 
(V,p) in the 2D region (x.c) = [0, 1][0, 1]. A 3D ver- 
sion of the problem corresponds to Koseff & Street ex- 
periment [1] with the domain (x,y,z) = [0, 1][0,3][0, 1]. 
The boundary conditions are: V = («.v,w) = (1,0,0) 
at the driven lid {z = 1);' V = 0 on the other walls, and 
eq. (5) holds on all the walls. The undefined constant 
in the pressure field is eliminated by setting p - 0 at 

z) = (0,0,0). 

2) Backward facing step flow (2D). We solve for a 
flow in a plane channel with sudden (step) expansion 
that corresponds to the experiment by Kim et al. [14]. 
The ratio of a step height H to the channel outlet width 
L is H /L= 1/3. The Reynolds number for the exper- 
imental flow is Re = 4.4 ■ 10 4 with H as a length scale 
[14], or Re L = 1.32- 10 5 in [15] when the channel outlet 
width L is used as a length scale. The flow profile V(v) 
is given at the inlet and p — 0 and stress free conditions 
at the outlet. We have used the following domain geom- 
etry: H - 1 , L = 3, total length Lj — 12, and the length 
from the inlet to the step L* = 4. 

3) Thermal convection in a cylinder (3D). Ther- 
mal convection flow in a differentially heated horizon- 
tal cylinder is considered. The experimental data by 
Bogatirev et al. [16] have been used for comparison. 
These data have been obtained during the ground tests 
for this device, a thermal convection sensor designed 
for the space micro-acceleration measurements; it was 
later flown in space on Mir station. Numerical results 
from the 3D finite volume simulations by Bessonov [17] 
have been also used for a comparison. The temperature 
distribution on a cylinder wall was (i) linear tempera- 
ture profile, and (ii) computed using a real, finite wail 
conductivity (adjoint problem). The body torce in (1) is 
F = (0.0 ,/toPrO), and Re~ x in (1) is to be replaced by 
Prandtl number Pr — v/k. Experimental data are avail- 

from 10^ to 


the dimension of a product xv is a square of length. We 
introduce a regularization length scale /, / 2 = xv and 
rewrite x“ as r - ^L^Rc = K • Re , where K = l 2 /L 2 . able for the Rayleigh number Ra in a range 
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10 5 that present the temperature difference between two 
points inside a cylinder. 

5 Finite element model and solu- 
tion method 

Proposed model allows efficient implementation by the 
finite element method. Both velocity and pressure arc 
approximated by the same order finite element basis 
functions. We used linear basis functions on triangle fi- 
nite elements in 2D and trilinear on hexaedral elements 
in 3D. 

Note that when t" = 0, we obtain the NS equations, 
and the FEM equations become incompatible. It is a 
well known problem, and can be resolved if a different 
order FE basis is used for the velocity and pressure ap- 
proximation (see e.g. [18]). When x* is in the range of 
10“ 7 to 10 -4 , our solution coincides with the NS one 
up to Re = 10 4 in a case of driven cavity problem. The 
value of r smaller than 1(T 8 results in the wiggles in 
the velocity profiles. At low Re , Re < 1000 the solu- 
tion practically coincides with the NS one for any T 
between 10“" 7 to 10” 2 . 

A CFD code FEMINA/3D [18], [19] modified for a 
new model was used in the numerical experiments. 
Three scalar equations (or four in 3D case) from (1), 
(4) are solved simultaneously. The nonlinear algebraic 
system of equations, resulting from the FEM discretiza- 
tion. is solved iteratively by a Newton method. Cor- 
responding linear system is solved by robust precondi- 
tioned iterative CGS-type method with preconditioning 
by the high order incomplete decomposition. We used 
a CNSPACK linear solver with a compact sparse matrix 
storase (see on details of the solver in [20],[22]). A first 
order decomposition was found to be quite enough tor 
the RNS, while the second or third order is required for 
the NS related linear systems [20]. 

The solution is considered converged when an alge- 
braic equations residual norm is || r ||oo< 10 4 (typi- 

cal initial residual value were 1 to 10"). In the case 
of the thermal convection problem, the energy equation 
is solved separately at each time step or Newton iter- 
ation. Selected solution method and a corresponding 
software implementation make it possible to soh e large 
time-dependent problems (up to 300,000 unknowns) on 
the SGI 02 workstation. 

6 Results and discussion 

Initially for each problem we varied the value of x or K 
and compared numerical results with the experimental 


ones and with the NS. At low Re < 1000 the solution 
practically coincides with the NS one for any x’ be- 
tween 10 -7 to 10 -2 for problems 1 and 2. We solved 
the transient problem equations (1), (4), or (1), (4) and 
(3) and have been able to obtain stationary solutions. A 
more efficient approach in such a case was a solution 
of the stationary problem by the Newton technique. In 
those cases there is a good agreement between the NS 
and the experiment when it is available (for example, 

Re = 50 to 500 in [21]). 

Mass conservation in the continuity equation was 
thoroughly analyzed. The local values of div\ of the 
numerical solution were examined. Note, that the nu- 
merical solution satisfies the discretized equation ex- 
actly only at specific points (nodes, integration points or 
collocation points etc. depending on the method used), 
and only approximately at an arbitrary point of the do- 
main. As the mesh size decreases, the accuracy of the 
approximation of the equation improves (typically as 
A 2 , the square of the mesh size). 

The value of divV for the numerical solution at an 
arbitrary point is non-zero, and we found that it was of 
the same order both for the RNS and NS numerical solu- 
tions. The values of divV exhibited spatial oscillations 
with amplitudes of 10~ 3 to 10" 2 . The averaged value 
of divV (integral over the region near the vicinity of any 
point) was the same order for RNS solutions as for the 
NS ones. The integral over the whole domain of divV is 
of the same order for both RNS and NS. Therefore, we 
see that in the numerical solution, a regularization term 
in the RNS continuity equation does not create more 
numerical mass fluxes, than in the case of NS. 

For a comparison we also present the results ob- 
tained with a standard Jt — e turbulence model, using 
the commercial code CFD2000 [4], A k - e turbu- 
lence model was developed for computations of tur- 
bulent flows, when the NS equations can not provide 
reasonable results. Viscous term in the NS momentum 
equation is replaced by a different stress term with a tur- 
bulent Reynolds number Re , . Two more nonlinear equa- 
tions are added to the system, that govern the turbulent 
kinetic energy k, dissipation rate e and provide evalua- 
tion of the Re,. These equations contain eight empirical 
constants for a standard model ( 1 1 or more in improved 
models). The standard k—t turbulence model is im- 
plemented in many commercial CFD codes and widely 
used in the solution of engineering problems. 

6.1 Driven cavity problem 

We used 81x81 nodes with a homogeneous triangular 
mesh for 2D computations. Same mesh with quadratic 
(6-noded) triangles and linear basis FE functions for 
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Figure 3: Driven cavity problem. Re = 3200. Compar- 
ison of horizontal velocity profiles (1-4) tor numerical 
(solid and dashed lines) and experimental (squares) re- 
sults and vertical velocity profiles (5-8) for numerical 
(solid&dashed lines) and experimental (triangles) re- 
sults : 1-NS, 2 - k — £ model, 3-RNS (2D),4-RNS (3D); 
5-RNS (3D), 6-RNS (2D), 7 - k - e model, 8-NS 



Figure 4: Driven cavity problem. Re = 10.000. Com- 
parison of horizontal velocity profiles (1-3) tor nu- 
merical (dashed lines) and experimental (squares) re- 
sults and vertical velocity profiles (4-6) for numerical 
(dashed lines) and experimental (triangles) results: 1- 
NS, 2 - * - e model, 3-RNS (2D); 4-RNS (2D), 5 - k - e 
model, 6-NS 


pressure and quadratic tor velocity was used for the NS 
computations. We note at once that further mesh refine- 
ment did not influence the results tor both flow regimes 
considered. For Re > 10 3 the RNS and NS solutions 
become noticeably different for x’ > l O'* 2 . We varied 
the value of K to match the experimental velocity pro- 
files at Re = 3200 presented in [1]. Good agreement 
was obtained, when K zz 1.5 ■ 10“-. 

Using the same value of K as above, we computed the 
RNS solution at Re = 10 4 , for that [1] also provides the 
data. The results are presented in Fig. 3. The error es- 
timates for experimental points are 1 to 10%, according 
to [1], Therefore, we can conclude, that the RNS so- 
lution is in better agreement with the experimental data 
than the NS solution. Indeed, at this value of Re there is 
sienificantly greater difference between the NS solution 
and the experimental results. The k — £ model results 
also disagree with the experimental data. 

Fis. 3 and 4 present the experimental measurements 
and numerical results for u and u- velocity components 
obtained from the RNS (2D and 3D), 2D NS, and 2D 
fc — £ model. The experimental data are shown for the 
symmetry plane (x.z) at y= 1 .5 along the lines x =0.5 
and z =0.5 of the 3D cavity, that has relative (x:y:z) 
dimensions 1:3:1. Since strong instantaneous velocity 
fluctuations were observed, the experimental points cor- 


respond to time averaged values of the measured veloc- 
ities [1]. 

The results obtained with RNS model are in good 
agreement with the experimental data for Re=3200 and 
10,000 and are an improvement upon previous results 
obtained using the Navier-Stokes equation. 

An important feature of our numerical experiments is 
the parameter K = / 2 /L 2 . We can estimate /, using the 
value of K found and the experiment description from 
[1]: L= 15cm, so / ~ 0.58 mm. This value of / is a good 
approximation of the "Kolmogorov microscale" l exp ~ 
0.5mm (the smallest scale of the flow nonhomogene- 
outy, viscous dissipation scale) that was observed in the 
experiment ([1], p.398). 

That is one of the ways to determine the actual value 
of r or K for the RNS model in advance, that is the 
experimental approach. We remind the reader that for 
low Re numbers the value of T* is actually insignificant 
and can be chosen in a rather wide range. It is the case 
for the problem 3 (thermal convection). 

For 3D flows we have solved the RNS for Re = 3200 
with the mesh of 41x41x33 nodes (221,892 unknowns), 
refined near the walls, for a half of the cavity. The sym- 
metry condition was used on a symmetry plane. Our 
results of modeling are presented in Fig. 3 as well. One 
can see that velocities obtained for both 2D and 3D pro- 
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Re=44,000 VELOCITY PROFILE 



Figure 5: Backward facing step flow, Re — 44,000. 
Comparison of the horizontal velocity profiles at x = 
5.33 for the experimental [14], RNS and k — z model 
results. 

files are close to the experimental data except in one re- 
gion. For the Re = 1 0 4 we did not see any 3D numerical 
results published for this case. The numerical solution 
for our model does not exhibit strong oscillations as a 
measured experimental velocities for both flow regimes 
considered. In the case of Re = 3200 we are even able 
to obtain a stationary solution by the Newton technique. 

6.2 Backward facing step flow 

The numerical results for a 2D backward facing step 
flow at Re = 4.4 ■ 10 4 (referenced as Re = 1.32- 10" 
in [15]) were obtained and compared with experimen- 
tal measurements [14]. We used a 1 10x60 mesh re- 
fined near the walls and did computations starting at 
low Reynolds number and increasing Re by small incre- 
ments until 44.000. Fig. 5 presents the computed ve- 
locity profile at x = 5.33 (jc = 0 at the begin of the step), 
the experimental mean velocity measurements [14] and 
the k — E model results. 

The RNS model results agree with experiment for 
both the velocity profile and the recirculation zone 
length X r . Our result is X r fH = 7.50 and the experi- 
mentally obtained one is X* xp /H « 7 ±0.5. The value 
of x' used in the computations was in a range from 10 
to 10“'. It did not influence noticeably the results. For 
a smaller r is was more difficult to reach the steady 
state solution at Re = 44.000; we have to use more and 
smaller increments in Re number to reach the final value 
of Re. If the increment in Re number was large, the flow 
pattern bifurcated to the unsteady flow with the vortexes 
periodically originating from the recirculation zone and 
flowing downstream. Our conclusion is that the numer- 
ical RNS solution for this problem does not depend on 


the value of T\ 

The solution with a standard h-Z model gives the 
velocity profile at x — 5.33 that has no backward flow 
at all [15]. A standard k-z model underpredicts the 
reattachement point (recirculation zone length X r ) by a 
substantial amount of 20-25% according to paper [15], 
where different turbulence models have been analyzed 
for this problem. 

6.3 thermal convection in a differentially 
heated horizontal cylinder 

The RNS model is applied to the solution of the problem 
of convective 3D flow in a differentially heated horizon- 
tal cylinder. Initially a linear temperature profile was 
assumed to be given on a cylinder wall. Experimental 
data by Bogatirev et al. [16] and finite volume simu- 
lations by Bessonov [17] are used for the comparison 
at Rayleigh number in the range 10 3 to 1.2* 10 s (Fig. 
6a). For the linear temperature distribution on a cylin- 
der wail we obtained a good agreement with the numer- 
ical results by Bessonov [17], obtained by the 3D finite 
volume method for the NS equations. The agreement 
with the experimental data was not good for Ra number 
more than 2000 (solid line in Fig. 6b). The computa- 
tions have been also done in [17] for a real, finite wail 
conductivity. A thermodiffusivity data for a stainless 
steel was used and adjoint problem was solved. The 
agreement with the experimental data has been signifi- 
cantly improved (dashed line in Fig. 6b). 

The value of t’ used was 10~ 7 to 10“\ and it did not 
change noticeably the RNS results. 

Another successful application of the RNS model to 
a 3D thermal convection flows with vibrations applied 
is presented in [23]. The application to the problem of 
magnetic field suppression of the semiconductor melt 
flow in crystal growth is considered in [11], where a 
good agreement with the theoretical asymptotic solution 
is obtained. 

7 Conclusions 

A numerical model for incompressible viscous flow is 
proposed. It is based on the regularization of the Navier- 
Stokes equations. Good agreement with the experiment 
is obtained for a driven cavity flow at Re - 3200 and 
10 4 . The regularization parameter found to be related 
to the experimentally observed spatial fluctuation scale. 

The RNS model solution approaches the NS one for 
small Re number. For high Re number the RNS solu- 
tion is more close to the experimental data than the NS 
solution (driven cavity flow, backward facing step flow). 
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Good agreement with the experimental data for a 3D 
thermal convection in a horizontal cylinder has been ob- 
tained at Ra = 10'' to lO 5 . 

The optimal value of the regularization parameter t" 
still remains an open problem. We have found that the 

I olu numerical solution depends on the regularization length 
§*§ l for the driven cavity problem, and the optimal value 
Smo of l is the same for both the 2D and 3D RNS solutions 
TEMP and for a11 the floW re S imes considered (Re = 3200 and 
r 10,000).’ The value of / agrees with the ’’Kolmogorov 
microscale” observed in the experiment [1]. 

The value of x* (or /) is not significant for a 3D 
thermal convection in a horizontal cylinder and for a 
backward facing step flow ( Re — 44,000). The RNS 
solutions agree with the experimental data and do not 
change even if t* is changed by few orders of magni- 
tude. 



a Experimental Data 
* R N S 


0 2 


6 8 10 
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Abstract 

A hierarchical family of methods for highly accurate solution for viscous flows at high Reynolds number, flows 
with boundary layer is presented. Note, that thin boundary layer can occur even at low Re, e.g. in flows of electrically 
conductive fluids under strong magnetic field at Hartman number Ha » 1 . The main feature of the methods proposed 
is the restoring of conservation laws at the boundary. Numerical solutions using this new approach compare favorably 
with available exact or asymptotic ones and experimental data. 


Introduction 

Difficulties with the numerical solutions of Navier-Stokes (NS) equations for high Re number flows have been ref- 
erenced usually to insufficient mesh resolution, complicated flow physics, turbulence etc. Shishkin et al. in senes 
of papers (1995-1997) shows theoretically that grid methods perform poorly in dealing with the boundary layer and 
provides the estimation for numerical solution error as 0(1) for uniform meshes [1]. It was proposed the construction 
of special meshes as a remedy for these difficulties with more optimistic error estimation as 0{h ) h is a mesh 

size with m = 9 or more. Failures of numerical solutions to obtain good agreement with expenmental data during 
last decades can be referenced as confirmation of Shishkin estimates (lid-driven cavity flow: expenments [2] versus 
simulations [3]; thermal convection in a vertical slot: experiments [4] versus [5], etc.). Successful results with imple- 
menting Shishkin-type strategy are obtained in [6], To construct Shishkin mesh, one has to use the boundary layer 
thickness that is not known in advance. Meshless methods have similar drawbacks in such problems (see e.g. 17]). 

The facts above and recent experiments with different formulations [8] and regularizations [10) for the viscous 
flow problem resulted in uncovering the main reason for poor numerical results: conservation laws are violated at the 
boundary in numerical treatments. The violation is done explicitly in FEM by dropping the weak formulation at the 
Diriehlet boundary (there is no choice); the same is not explicit but present in traditional FDM, spectral or multiqua nc 
(MQ) collocation methods (governing equations are not used at the boundary). That can be ignored tor we "P° se 
problem, but becomes crucial for a singularly perturbed boundary value problem, e.g. viscous flow with boundary 

^^We present the new alternative numerical solution methods that are more accurate, efficient and less mesh depen- 
dent These methods, raneins in the increasing of accuracy (and complexity of implementation), restore the conserva- 
tion laws at the boundary^ (i) a regularization approach, the RNS [10], (ii) a projection approach, (m) an extension of 
basis function space (for FEM or MQ[15]). 

Higher accuracy numerical solution methods for boundary layer 

A regularization approach to solving the Navier-Stokes equations 

We consider a flow of incompressible viscous fluid in a closed 2D or 3D domain Q. The Navier-Stokes momentum 
and continuity equations are: 

— + i-Vp - F = 0 


( 1 ) 



( 2 ) 


V-V = 0 

where Re = V 0 L/v is the Reynolds number, V 0 is the velocity scale, L is the hydrodynamic length sealery is the 
kinematic viscosity, and F is a body force. In the case of thermal convection the body force is F = GrRe •© <* 
(Boussinesq approximation), where © is a nondimensional temperature, Gr is the Grashoff number and e, is the unit 
vector in the direction of gravity. The energy equation is: 


39 


+ (vv)e = /v _l fl<r'v 2 e 


3t 


(3) 


where 9 represents nondimensional temperature, scaled by 9 (I - T co ij)/AT with AT = 7),,,, -T t n ij. The Prandtl, 
Grashoff and Rayleigh numbers are respectively Pr = v/k, Gr = Ra/Pr, and Ra - $ATgL k v , where [i, g, are 
the coefficients of thermal expansion, gravitational acceleration, and of thermodiffusivity. 

Boundary and initial conditions complete the problem formulation and they will be given for each pro em. 

A proposed regularization consists in modifying the continuity equation (2) to become 


V • V = T*V • (Vp — F) 


(4) 


► 0 eq. (4) approaches the continuity equation (2). The boundary 


where X* is a small regularization parameter. For x* 
condition for the pressure on a wall is 

(Vp — F) n = 0, (5) 

where n is a unit vector normal to the wall. Equations (4) and (5) present the basis of this method. Equation (5) ensures 
a balance of the components of the forces that are normal to the region boundary. 

First results with a more complicated model have been presented in [8] and [9]. Further numerical experiments 
have shown that accurate results can be obtained with the simpler model [10]. In [10] we have shown that equations 
(1), (3), and (4) with (5), called as the regularized Navier-Stokes equations (RNS model), give a better agreement with 
the experimental measurements for high Reynolds number flows than the traditional NS solution by FEM. 

The pressure Laplacian term on the right-hand side of eq. (4), has also been obtained by Lohner as a result ot 
consistent treatment of the time-advancement for the continuity eq. (2) [11]. There a discretization time step Ar is 

used in place of our parameter x* and the boundary condition for the pressure was not specified. . 

Tire regularization parameter f is expressed dimensionally as x* = tL~%, where t is time The dimension of w 
is leneth squared. We introduce the regularization length scale/ with l 2 = tv and rewrite x asx =IL Re - K ■ Re. 
where K = I 2 !L r. The optimal value of X (or /) is not known in advance. We found that for problem with smooth 
boundary condition the solution is undependable of the value of x*. with x* is the range from 10 to l . 

Such a regularization has an additional useful feature for the FEM. It allows to use the same order finite element 

approximation for the velocity and pressure . . nn . u t - i 

Numerical results that we obtained for 2D/3D flow problems are in dramatically better agreement with analytical 
solutions and experimental flow measurements. The numerical 2D/3D solutions by FEM with this strategy employed. 

called RNS against the 3 D driven cavity flow data [2] and results by other methods are shown in Fig. . e resi ua s 

of numerical solutions, shown in Fig. 2(a), demonstrate by one to two order residual reduction in the boundary layer 
for the Navier-Stokes momentum equation. We recognized that eq. (4) is an approximation of the eq. ( ) in projec ion 

Further development of the method can be the use of the governing equations at the boundary. This task is 
complicated in the frame of the low-order FEM scheme; a projection approach can be implemented wi^ Hermne 
finite elements (not presented here). This can be done simpler in the frame of the multiquadnc method (MQ) [ 12]. 

Another approach, that is being developed for the FDM, involves high order finite difference schemes that are 
specially constructed to deal with the boundary layer solution [16], [17]. 

PDE collocation at the boundary 

To have the PDE to be satisfied at the boundary- one needs to extend the basis function space in the MQ by N h functions 
for each PDE. where N„ is the number of boundary nodes. It can be done by adding the layer of nodes adjacent to the 
boundary and related MQ basis functions, thus allowing to add a set of equations obtained via collocation of the PDE 

0,1 SS ID singular perturbation problem from [7] by multiquadric global basis function method [12], using 
this strategy with 21 nodes, versus exact solution is shown in Fig.2(b) for e = 10 ;41 nodes are enough lor e - 10 


2 



(not shown), 
method with 


The gain in the accuracy with this strategy is two to three orders compared to the original Kansa-MQ 
101 nodes[151. 
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(b) backward step flow, Re = 44,000(132,000), V x X 5.33 


(a) Cavity, Re = 3200, V x at X = 0.5 


Figure 1: (a) Lid-driven cavity problem. Re = 3200. Comparison of horizontal velocity profiles (lines 1 to .4; I for 

numerical solution and experimental data (squares) by Koseff and Street (1984). 1-NS. 2 - ■ emo e , . ’ 

4-RNS (3D); (b) Flow over a backward facing step. Re = 44,000(132,000). Comparison o t e onzonta ve oci \ 
profiles for the numerical solutions with RNS [10], k - e model [14] and experimental data (squares) [13] at x - 5.33. 

Summary 

Proposed higher accuracy solution methods are used in applications to 2D and 3D flows at high Reynolds number 
and flows with thin boundary layer, and compared favorably with experimental data and asymptotic solutions (when 
available): flow in channels, thermal and vibrational convection (Fig. 3, [18]), electrically conductive fluid flows un er 
strong magnetic fields etc. [10]. 


LID-DRIVEN CAVITY FLOW Re-3,200 


EXACT AND NUMERICAL SOLUTIONS 
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(a) residuals of NS and RNS solutions 


o 0 02 0 04 0 06 0.08 

COORDINATE X 

(b) ID boundary layer solution 


Fieure 2: (a) Residuals of numerical solution in the X-momentum equation at c 0.5 for lid-driven cavity problem. 
p/ = 3 200 for the 2D NS solution and RNS solution on the same 81*81 mesh. Residuals are estimated by differ- 
entiation of numerical solutions, using interpolating polynomials. The RNS residual is lower by one to two order 
than the NS one; (b) ID problem with boundary layer: eh„ + k, = j>, «(0) - 0 n(l) - 1, comparison o "umenca 
(multiquadric method. 21 nodes, strategy (iii) and exact w = (1 -e x/ )/(l -e ) solutions for e 10 . 

in [0,0.1], 
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Figure 3: Typical instantaneous 3D thermovibrational melt flow patterns for: (a) translational vibration at Og, Ra 0, 
p av _ 7 25 . io 4 , Pr = 0.01, co = 1 00//'. lateral vibration is applied along x-direction (horizontal), velocity compo- 
nents Vx Vx V- pressure P. temperature T, temperature disturbance AT, and velocity module; 

(b) circular' polarized vibration, Ra = 7.25- 10 3 , Ra v = 7.4- 10 6 , Pr= 15. to = 10 Hz : velocity components 
temperature T, concentration C, velocity module, and temperature disturbance AT. 
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Abstract 

We present an efficient solution technique for large sparse nonsymmetric algebraic linear 
systems, related to the coupled solution of incompressible viscous fluid flow equations. An 
iterative solution technique with high order incomplete decomposition as preconditioning is used 
in this method. A developed CNSPACK linear solver for flow problems on 2D/3D unstructured 
meshes is briefly described. Analysis of the efficiency of the proposed approach is demonstrated 
on 2D and 3D flow problems. Numerical experiments show that the computational complexity of 
the proposed method appears to be 0(N 5/4 ). 

Introduction 

Numerical solution of incompressible viscous flow problems by the finite-element method leads 
to lar°e linear systems of equations with sparse non-symmetric matrices. For the Navier-Stokes 
equations in primitive variables with the continuity equation, treated directly, the system of 
related linear equations combines many unwanted properties of algebraic systems. It is non- 
svmmetric. not positive definite and has zero diagonals in rows related to the discretized 
continuin' equation. The advantage of the coupled solution is in the improved stability of the 
numerical solution for stationary flow at high Reynolds numbers. Also, much larger discrete time 
steps for transient flow analysis can be used. The drawback of the coupled solution is that it is 
expensive in terms of computer resources, such as memory and computation time. 

Recent advances in iterative solution techniques resulted in more efficient solvers, based 
on CGS and BiCGS [1.2] and GMRES methods [3,2]. A preconditioning technique can be used 
to accelerate the convergence of iterative solution, thus resulting in dramatically fewer iterations, 
tvpicallv 10 to 100 instead of the order of 1000 iterations without preconditioning. 

' 1 Recurrent formulas have been published to construct a preconditioner for structured (3- or 
5-dia2onal) matrices [4], Simple first order preconditioners are typically used for unstructured 
matrices that are related to unstructured meshes for complicated domain geometry. The efficiency 
of such preconditioning is low and, in a case of the incompressible Navier-Stokes equations, 
solved with quadratic for velocity and linear for pressure finite elements, such a preconditioning 
mav result in a divergence of iterations. In [5] we proposed the high order preconditioners for 
unstructured matrices, which improve convergence of the iterative solution process. In this paper. 



we analyze the efficiency of this approach on 2D and 3D flow problems. NumericaUxperiments 
show that the computational complexity of the proposed method appears to be 0(N ). 

Problem formulation 


The following are governing equations for viscous flow in closed domain G. The momentum 
equation is 

— + (UV)U- — AU + Vp = 0 (!) 

dt Re 

while the continuity equation is 

div{ U) = 0. W 

Here U is velocity vector, p is pressure, and Re is the Reynolds number. We present the following 
examples: (i) flow in a square cavity with driven lid; (ii) viscous flow around the system of 
circular cylinders; (iii) a flow in ion-exchanger channel (Fig. 1), and (iv) flow in a 3D cavity (Fig. 

2 ). 



Fie. 1 : Test problems, meshes and flow patterns: (a) driven cavity, 2D, (b) flow around the system 
of circular cylinders, and (c) flow in ion-exchanger channel. 

Boundary conditions are: (i) velocity is equal to zero at the bottom and side walls, and a 
horizontal component a of velocity is equal to 1 on a driven top lid, (ii,iii) inlet ve ocity pro i e 
u—u( v j, V =0 is given at the left domain side, zero velocity on top and bottom, and at cylinder 
surfaces, and p = 0 and “soft” boundary conditions for velocity at the channel exit,at the rig t 
side of the domain. Initial conditions are zero fields for all variables. 

Solution method 

We used the finite element method on 6-noded triangles with quadratic for velocity and linear for 
pressure finite elements for 2D problems. A new regularization method for viscous flows at high 



Reynolds numbers with trilinear interpolation for both velocity and pressure was used in the 3D 
driven cavity example [7,8] with a geometry corresponding to the experiments [16]. 

The discretized linear algebraic systems are solved using the CGS-type iterative method 
with preconditioning by the incomplete decomposition of the original matrix. Comparing the 
CGS and GMRES methods [1, 3] in our examples, we found that both methods convergence well, 
if a good preconditioner is used. The CGS method needs comparably less memory to store only 

eight work vectors. •' 

To reduce the memory requirements, a compact storage scheme for matrices was used. 

We stored only the nonzero entries. The incomplete decomposition (ID) was a product of 
triangular and diagonal matrices, P = LDU. To avoid a diagonal pivot degeneration we use the 
Kershaw diagonal modification [6]. If the value ot diagonal element was small, i.e. 

I a ..|< a = , the diagonal was replaced by a . Here cr and /iare the maximum magnitudes of 

row and column elements, and 2* is a machine precision (f bits in mantissa, see details in [6]). 



Fia.2: 3D driven cavity problem: mesh and horizontal velocity component is shown (dark gray 
color shows larger velocity component magnitude). 

For the first order ID, the matrix P has the same non-zero entry pattern as the original matrix. For 
a second order or higher ID. matrix P has one or more additional entries near the locations of the 
non-zero matrix entries, where the original matrix entries are zeros. 


Numerical experiments 

The computations were done on the SGI 02 machine for different mesh resolutions and different 
incomplete decomposition orders. Iteration termination criteria was the value of the module ol the 
residual: r < 10' 6 . 10' 8 , or 10 12 . The solution times for a driven cavity problem were compared 
for different solution methods (Fig. 3). 

The results of our numerical experiments can be summarized as follows: 


(1) The iterative solution without a preconditioning, ID = 0, is time-expensive, the number of 
iterations is of the order of the number of unknowns. The memory requirements are smallest. 

(2) A simple incomplete decomposition of the first order, ID = 1, can result in the divergence of 
the iterations. That may be a reason, why this approach is not widely used yet. 

(3) The incomplete decomposition of the second or third order can reduce the number of 
iterations dramatically, e.g. for N up to 10,000 the number of iteration did not exceed 200, 
and, for most cases, just 50 iterations were enough for convergence. The memory 
requirements are about twice as high as for ID = 0. 



Fig. 3; Solution time versus number of unknowns N for 2D and 3D driven cavity problem, 
comparison of different software packages. FRONTAL - direct solution by frontal method. 
CNSPACK-2D- solution with proposed method, P2-P1 finite elements, matrices stored in 
double precision, incomplete decomposition of the second order; Y12M - solution with 
Y12M sparse matrix solver package; CNSPACK-3D (lower curve) - solution of 3 
problem; matrices are stored in single precision, a first order incomplete decomposition is 
used. Theoretical asymptotics for direct solution and proposed method, the N , and N 
lines, are plotted for comparison. 

(4) A new regularization approach for viscous flows, that works well with thin boundary layers 
[7.8], was used in our tests for 3D cavity problem. It resulted in a much faster solution 
convergence. For N up to 300,000 the number of iterations to converge was of the order of 30 
to 50. And in this case, the first order incomplete decomposition worked quite well. 

(5) We compared the solution time for the same linear system by different methods: direct 
solution by frontal method [5], iterative solution with Y12M sparse matrix package [15] and 
by our CNSPACK software (Fig. 3). We have been unable to use the SPARSKIT [1 1] and 
SPARSE [12] solvers, as those failed for zero matrix diagonal elements. 




(6) Convergence rate of the proposed method does not reduce with increasing the Reynolds 
number Re (from 10 to 1000). We observed even faster convergence for larger Re. 

(7) Proposed method can solve the linear system with a high accuracy, much higher than obtained 

by a direct solution. w 

(8) The node/equation numbering is of great importance for the rate of convergence [13, J. e 

found that the optimal node or equation renumbering is quite necessary for unstructured 
meshes, and ID = 2 should be used for fast convergence. 5/4 

(9) It appears that the computational complexity of the proposed method is 0(N ) (see Fig. 3). 

Aspects of the efficiency of the method on superscalar pipelined microprocessors and parallel 
architectures were considered. We found the solution time can be reduced by a factor of three by 
using special optimization techniques for a superscalar pipelined architecture. 

Proposed method is used in FEMINA/3D code that is used for solving 2D/3D incompressible 
viscous flows, 3D thermal convection, magneto-hydrodynamics flows, and 3D thermal 
vibrational convection in Bridgman crystal growth configurations [7-10], If the number of 
unknowns is in the range of N- 50,000 to 500,000, the problems are solved on a low-end SGI 
workstation or PC, using our FEMINA/3D code (for N of the order of 10 , the supercomputer is 
needed). For example, we solve large 3D time-dependent problems of thermo-vibrational 
convection with 300.000 unknowns on the low-end SGI 02 workstation within a couple of hours 

[8],[9]. 

Conclusions 

We have presented an efficient solution technique for large sparse nonsymmetric algebraic linear 
systems for coupled solution of incompressible viscous flow equations. An iterative solution 
technique with high order incomplete decomposition as preconditioning is used in this method. 
Analysis of the efficiency of this proposed approach is demonstrated on 2D and 3D flow 
problems. Numerical experiments show that the computational complexity of the proposed 
method appears to be 0(N 5/4 ). 
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Abstract 

The RNS method was used to solve viscous flow problems at high Reynolds numbers with a thin boundary layer. 

This method is related to the generalized hydrodynamic equations proposed by Alexeev (1994). It can be interpreted 
as a regularization of the Navier-Stokes equations. Numerical solutions using this approach compare favorably with 
experimental data. The method is discussed and numerical solutions are compared with the experimental data for 
a 3D driven cavity flow at Re = 3200 and 10,000, 2D backward facing step flow at Re = 44,000, 2D channel flow 
at Re number up to 10 6 , and a 3D thermal convection in a cylinder at Ra = 1000 to 150,000. Comparison with the 
analytical asymptotic solution is provided for a thermal convection, in the electrically conducting fluid suppressed by 
a strong magnetic field at Hartman numbers Ha up to 20,000. This proposed model is not a turbulence model, and no 
additional equations are introduced. 

Keywords: High Reynolds number flows, Alexeev equations, Regularization, Finite element method. 

1 Introduction 

Difficulties with the numerical solutions of Navier-Stokes (NS) equations for high Re number flows have usually been 
referenced to insufficient mesh resolution, complicated flow physics, turbulence, etc. In series of papers, Shishkin et al. 

( 1 995- 1 997, see e.g. [ 1 ]) demonstrated theoretically that grid methods perform poorly when dealing with the boundary 
layer and provided an estimation for numerical solution error as 0(1) for uniform meshes [1]. Asa remedy for these 
difficulties, the construction of special meshes was proposed, with more optimistic error estimation, as 0(h ), h is 

the mesh size, with m = 7 or more. During last decades the failure of numerical solutions to obtain good agreement 
with experimental data confirms Shishkin estimates (lid-driven cavity flow: experiments [2] versus simulations [3]; 
thermal convection in a vertical slot: experiments [4] versus simulations [5], etc.). Successful results implementing 
Shishkin-type strategy are obtained in [6]. The boundary layer thickness, which is not known in advance, is used, to 
construct Shishkin mesh. Meshless methods have similar drawbacks [7]. 

We have developed an alternate approach, that is more accurate, less mesh-dependent and can be interpreted as a 
regularization of the Navier-Stokes equations. It is based on the mathematical model associated with the generalized 
hydrodynamic equations by Alexeev [8], [9]. The numerical solutions agree well with the experimental measurements 
for a set of flow problems at high Reynolds number and for flows with thin boundary layer [10], [11], [12], [13], 

[ 1 4] [ 1 5]. This is not a turbulence model, and no additional equations are introduced. 

The model was successfully compared to: ( 1 ) 3D driven cavity flow data by Koseff and Street ( 1 982) at Re — 3200 
and Re = 10,000: (2) 2D backward facing step flow by Kim et al. (1980) at Re = 44,000; (3) 3D thermal convection 
in a cylinder at Ra = 1000 to 150,000; and (4) asymptotic solution for a thermal convection in the semiconductor melt 
suppressed by the magnetic field at Hartman numbers Ha up to Ha = 20,000. Numerical results for 2D channel flow 
at Re number up to 10 6 are also presented. 

2 A regularization approach to solving Navier-Stokes equations 

Governing equations. To consider a flow of incompressible viscous fluid, in a closed 2D or 3D domain ft, Navier- 
Stokes momentum and continuity equations are: 

^X + ( V V)V-/?e-'V 2 V + Vp - F = 0 (1) 

dt 
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(a) Re=3200. velocity U{:.) at x = 0.5 


(b) residuals of NS and RNS solutions 


(c) Re= 1 0.000, velocity U(z) at x = 0.5 


Figure 1' Lid-driven cavity problem. Re = 3200: (a) horizontal velocity profiles (1-4) for numerical solution (solid 
and dashed lines): l -NS, 2 - standard k - e model. 3-RNS (2D), 4-RNS (3D); squares - experimental data by Koseff 
and Street (1982) in the symmetry plane Y = 1.5. (b) residuals of numerical solution in the momentum equation (1) 
for the "’D NS and RNS solutions at .v = 0.5. Residuals are estimated by differentiation of numerical solutions, using 
interpolating polynomials, (c) Re = 10 4 . numerical horizontal velocity profiles ( 1 -3) and experimental data (triangles). 
Ak-Z model solution is obtained with commercial code CFD2000. 

v-v = o 

where Re = V 0 L/v is the Reynolds number. V 0 is the velocity scale. L is the hydrodynamic length scale, y is the 
kinematic viscosity, and F is a body force. In the case of thermal convection the body force is F - GrRe -Q- e, 
(Boussinesq approximation), where 0 is a nondimensional temperature, Gr is the Grashotf number and e K is the unit 
vector in the direction of gravity. The energy equation is: 


^ + (VV)e = Pr~ l /?e~'V :! e (3) 

3t 

where 0 represents nondimensional temperature, scaled by 0 = {T - T a ,ij)/AT with AT = T hm ~ ^ Prandtl. 

Grashoff and Rayleigh numbers are respectively Pr = v/k, Gr = Ra/Pr , and Ra = pATgL-V v~>, where p, g, k are 
the coefficients of thermal expansion, gravitational acceleration, and of thermodtffusivity. 

We have analyzed the generalized hydrodynamic equation, proposed in [8] (a review in [9]), for the case of incom- 
pressible viscous flow and kept only a few main order terms, spatial fluctuations, in the continuity equation. This may 
be interpreted as a regularization of the Navier-Stokes equations. 

A proposed regularization involves modifying the continuity equation (2) to become 

V-V=x'V-(V/7- F) ( 4 ) 


where X 9 is a small regularization parameter. For T* 
condition for pressure on a wall is 


■ 0. eq. (4) approaches the continuity equation (2). The boundary 


(Vp — F) ■ n = 0, (5) 

where n is a unit vector normal to the wall. Equations (4) and (5) present the basis of our method, which takes into 
account only a few of many additional terms of the generalized hydrodynamic equations, called fluctuations (tempora 
and spatial), in [8], Equation (5) is a condition of absence of the fluctuations on walls, according to [8]. 

Preliminary results with a more complicated model are presented in [10] and [1 1]. Further numerical experiments 
have shown that satisfactory results can be obtained with the simpler model presented above [12]. Numerical solutions 
of eqs. (1). (3), and (4), called regularized Navier-Stokes equations (RNS), with the boundary condition (5), give a 
better agreement with the experimental measurements for high Reynolds number flows than the traditional solution of 

Navier-Stokes eqs. ( 1 ) and (2). (NS), with the finite element method (FEM). 

Note that numerical formulations, containing extra terms in the momentum and continuity equations, have been 
proposed in the frame of kinetically consistent numerical schemes, developed by Elizarova and Chetverushkin [ ], 

[18]. The justification of introducing extra terms into hydrodynamic equations is discussed in [19] from a physical 
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(a) Flow over a backward facing step. Re = 44,000 


(b) 2D channel flow, velocity U{z) at .v = 2.5 


Figure 2: (a) Flow over a backward facing step. Re = 44.000(1 32,000). comparison of the horizontal velocity for the 
RNS solution. k — £. model [25], and the experimental data (squares) [23] at x = 5.33; (b) 2D channel flow, horizontal 
velocity profiles t/(z) at .v = 2.5. computed with RNS Re = 5 • lO 3 (laminar and turbulent), 10 , 10 and 10 . 

kinetics viewpoint. The pressure Laplacian and other terms in the discretized continuity equation, have also been 
proposed by Lohner as a result of consistent treatment of the time-advancement for the continuity eq. (2) [16]. 

Further! eqs. (1), (3), and (4) with the boundary condition (5) are treated as mathematical model, having a control 
(regularization) parameter X* . This regularization parameter t* is expressed dimensionally as X* = XL V 0 , where x is 
time. The dimension of xv is a length squared. We introduce the regularization length scale / with /' - xv and rewrite 
X* as x* = l 2 L- z Re = K ■ Re , where K = I 2 /L 2 . The optimal value of x (or /) is not known in advance. We found that, 
for problems with smooth boundary conditions, the numerical solution only slightly depends on the value of x . with 
X* in the range of 10" 8 to KT 2 . This proposed regularization has an additional useful feature tor the FEM. It uses the 
same order finite element approximation for velocity and pressure for RNS. 


3 Numerical experiments with RNS 


In this section, we present the numerical results of a few flow problems, and compare these with experimental data. 
The FEM was used to discretize the governing equations. Algebraic equations for momentum and continuity were 
solved simultaneously using the iterative CNSPACK solver [20], with preconditioning by a high order incomplete 
decomposition. Iteration termination critetion was a convergence of relative residuals to 1 0 (or 10 ). 

Driven cavitv problem. Shown in Fig. 1 are the numerical solutions, with RNS employed t^l.fU], agains 
the and 3D driven cavity flow data at Re = 3.200 and 10,000 [2], and results by other methods. We used 81x81 
node uniform triangular mesh for the RNS. The same mesh, with quadratic for velocity and linear for pressure t finite 
elements was used in a standard FEM solution of eqs. (1.2). These results are labeled as NS in Fig. 1(a) and 1(c). 
The RNS model parameter, X* (or l 2 in l 2 L~ 2 \ was varied to match one of the experimental velocity profiles. A 
value l 2 L~ 2 w I 5 - I0 -5 resulted in good aereement for all the velocity profiles and both flow regimes, e 
and 10.000. To our surprise, the dimensional value of / ~ 0.58mm was a good approximation to the experimentally 

observed "Kolmogorov microscale” l exp ~ 0.5mm (see [2], p. 398). , O01 a(n 

For 3D flow we solved the RNS at Re = 3200 for half of the cavity with a mesh of 41x41x33 nodes (221 .89- 
unknowns). refined near the walls. The symmetry condition was used on the vertical symmetry plane y - 1.5. lhe 
same 2D value of K = 1 5 ■ 10 -5 was used in the 3D computations, the modeling results are presented in Fig. as we . 
One can see that the velocities obtained for both the 2D and 3D profiles are close to the experimental data, except in 
one resion. We did not obtain the 3D stationary solution at Re = 10 4 . The solution went unstat.onary at Re > 8.450. 
According to Baggett and Trefethen [22]. the stationary solution exists, but the basin ol attraction ot this solution may 

be extraordinarily narrow, having a width of 0{Re a ) for some a<-l. . , 

Mass conservation in the continuity equation was thoroughly analyzed. The local values of div\ for the numeric 
solution were examined. The conclusion was that a numerical solution for the RNS model has the same order ot error 
in the mass conservation as the NS solution with the same number of mesh nodes [ 1 2]. 





(a) Temperature 


(b) AT ~T\-Ti versus Ra 


Figure 3: (a) Thermal convection in a 3D differentially heated cylinder; (b) comparison of the temperature difference 
AT versus Ra with the experimental data [27] and finite-volume computations: squares - experimental data, solid line 
- numerical results for perfect wall conductivity [26], triangles - RNS results; dashed line - numerical results when 
the steel properties for the wall conductivity have been used (adjoint problem);here AT is a temperature difference 
between locations marked as 1 and 2 in (a). Both the NS and RNS solutions are nearly identical here. 

The residuals of numerical solutions in the momentum equation (1) are presented in Fig. l b. It shows one to two 
order reductions for the RNS solution residual in the boundary layer, compared to the NS solution residual. 

Flow over a backward facing step. The numerical results for a 2D flow over a backward facing step of height 
H, H = LI 3 (L is a channel height) at Re = 4.4 • 10 4 (or Re L = 1.32 • 10 5 in [24]) were obtained and compared with 
the experimental measurements of [23], We used a 1 1 0x60 mesh, refined near the walls, and started the computations 
at a low Revnolds number. We raised Re in small increments until reaching 44,000. Fig. 2(a) presents the computed 
velocity profile, at x = 5.33 H (x = 0 at the edge of the step), and the experimental mean velocity measurements [23]. 

The RNS model output satisfactorily agrees with the experimental data for both the velocity profile and the recircu- 
lation zone length X r . We computed X r /H = 7.50. while X? p /H « 7 ±0.5 was obtained experimentally. The value of 
r used in the compulations was in the range of 10" : to 1CT 4 . This did not noticeably influence the results. For smaller 
values of x\ it was more difficult to reach the steady state solution at Re = 44,000; more and smaller increments in 
Re had tobe used. If the increments in the Re number were large, the flow pattern bifurcated to an unsteady flow, wit 

vortices periodically originating from the recirculation zone and flowing downstream. 

The solution with a standard k - e model shows the velocity profile at* = 5.33 that has no backward flow. A 
standard k-z model underpredicts the recirculation zone length X r by a substantial amount, 20-25% according lo 
[24] where more sophisticated turbulence models have been proposed for this problem. 

Flow in a 2D channel, of height H= 1 and length L = 4, was the subject of a few experiments with the RNS at 
Reynolds numbers = 5 - 10\ 10 4 . lO 5 and 10 6 . An 81 x 100 mesh refined near the walls was used. Inlet flow profiles 
were (i) U = 1 and (ii) U = 6‘( 1 - z). We were obtained both parabolic and “turbulent” flow profiles for Re up to 10 . 
depending on the inlet flow conditions and the value oft*. To obtain a “turbulent” flow profile at Re = 5,000 with the 
inlet condition (ii), we started with t* = 0.01 . and were then able to keep this “turbulent ' profile type at reduce t . 
down to HT 4 (Fig. 2b). The boundary layer thickness is about 8~ Re~'i (obtained graphically from Fig. 2b). 

3D Thermal convection in a differentially heated horizontal cylinder. Alinear temperature profile is given 
assumed on a cylinder wall. Experimental data by Bogatirev et al. [27] and finite volume simulations by Bessonov 
[26] are used for comparison at Rayleigh numbers in the range of 10 to 1 .2 • 10' , Pr — 0.9 (Fig. 3a). or t e inear 
temperature distribution on a cylinder wall, we obtained good agreement with the numerical results [26], We used 
a 17.4357 node hexaedral mesh refined near the walls. Agreement with the experimental data was not good tor 
Ra > 2000 (solid line in Fig. 3b). Therefore computations were done for a real, finite wall conductivity [ ]. e 

thermodiffusivity data for stainless steel was used, and the adjoint problem was solved. Thereby, agreement with the 
experimental data has been significantly improved (Fig. 3b, dashed line). The value of X* used was 10 to 

which did not noticeably affect the RNS results. 

Magnetic field suppression of the semiconductor melt flow, modeling with RNS, is considered m [13] and [14], 
The application of magnetic fields is a promising approach for reducing convection during directional solidification 
of electricallv conductive melts. Current technology allows for experiments using very strong static fields, for which 
nearly convection free segregation is expected in melts exposed to stabilizing temperature gradients (vertical Bridgman 
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(a) Velocity V s (x), Ha = 2 • I0 4 


(b) ly zoom 
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Figure 4: (a)Vertical velocity profile V y (x) at y = 0.25 H, x in [ — 1, 1], Ha — 2.2- 10 4 ( B — 50 Tesla), (b) details 
of the velocity profile (a) in a boundary layer, .v from — 1 to —0.99; (c) Summary: maximum value of horizontal 
velocity magnitude versus magnetic field B for aligned and misaligned by 0.5 degree magnetic field and gravity 
vectors. Predicted theoretical asymptotic Umax ~ Ha~ 2 is observed starting from B = 0.5 Tesla (Ha = 220). 

method) [28]. The governing equations solved for this problem are eqs. (1), (3), and (4), where F. the body force due 
to the magnetic field (Lorenz force), is given by 

F = PrHa 2 [(\ x e B ) x e B \. 

Here F =(Pr(Ha)~U , 0) for the 2D case of a vertical magnetic field and U is the horizontal component of the velocity. 
The Hartmann number, given by Ha = LB 0 J%, is in the range of Ha = 100 to 10 4 , for the materials and magnetic 
fields under investigation. Here p.c are the density and electrical conductivity. Bo is the magnetic field intensity, and 

e B is the unit vector in the direction of the magnetic field. 

The computations are difficult because of the thin boundary layer, although the velocity of the generated flows is 
extremely low. Re ~ 10" 1 to lO -6 . A high value of the Hartmann number results in a relatively small coefficient at the 
highest derivative of the velocity in the momentum equation. Solutions of such problems exhibit thin boundary layers 
of thickness 5 ~ tfa" 1 , along with •‘equivalent” Reynolds number Re tqv ~ Ha 2 , Re eqv = 4 • 10 4 for B=0.5 Tesla, and 
Re . = 4 • 1 0 8 for B=50 Tesla. Some of the results for 2D models are shown in Fig. 4. The RNS numerical solution 
is rather smooth even for a very thin boundary layer, with thickness 8 ~ 10 -4 (Ha = 2 • 10 4 ). Other methods tested in 
[13], [14] (including industrial code) did not provide acceptable results or failed for Ha > 100. 


Conclusions 

The RNS method was review with applications to 2D and 3D flows at high Reynolds number and flows with thin 
boundary layers. The numerical results compared favorably with experimental data for driven cavity flow, flows in 
channels, thermal convection, and asymptotic solutions for electrically conductive fluid flows under strong magnetic 
fields The RNS method is used for modeling 3D thermo-vibrational convection in Bridgman melt configurations [15]. 
Similar ideas have been used successfully to improve the accuracy of the meshless multiquadric radial-basis function 
methods [29], [30]. 
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